Statistiques — Modèles linéaires

Notions fondamentales de la modélisation statistique

Population, échantillon

Les statistiques désignent un ensemble d’outils mathématiques issus de la théorie des probabilités permettant de décrire de façon synthétique des propriétés (des variables) d’un ensemble d’objets (un échantillon) tirés d’un ensemble plus grand et abstrait (la population). Ces outils visent en particulier à inférer des informations synthéthiques sur la population et à prédire les variables associées à un nouvel échantillon issu de cette population.

Une statistique est une quantitée calculée à partir des échantillons et supposée traduire une information sur la population. Par exemple, un histogramme.

Exemple fil rouge du cours - Un échantillon d’arbres dans la forêt de Grésigne en 2020, sur 9 placettes circulaires de 1 hectare au sein du massif. Plus d’informations sur le contexte scientifique d’acquisition de ces données sont disponibles à ce lien

dataIni <- read.csv("dataT3_simplif.csv")
dataIni$releve[dataIni$releve=="BLO_12b"] <- "BLO_12"
vLabRelFull <- names(
  sort(summary(as.factor(dataIni$releve)),decreasing=TRUE)[1:9])
dataIni <- dataIni[!is.na(match(dataIni$releve,vLabRelFull)),]

Décrire un échantillon

Visualisation de l’échantillon

Dans l’exemple fil rouge du cours, on va s’intéresser en particulier au diamètre (en cm) des arbres de la forêt de Grésigne. La variable diamètre est dénomée DBH (acronyme anglais pour Diameter at Breast Height).

N <- nrow(dataIni)
barplot(dataIni$DBH,xlab="Echantillons",ylab="DBH (cm)",names.arg=1:N,
        cex.names = 0.6,las=2)

Cette représentation est une première prise en main des données, mais elle ne permet pas d’appréhender la distribution de la variable d’intérêt au sein de l’échantillon. L’histogramme permet une visualisation de cette distribution :

hist(dataIni$DBH,main="DBH")

Attention c’est une représentation arbitraire, l’impression visuelle dépend des limites de catégories employées.

hist(dataIni$DBH,main="DBH",breaks=10:140)

hist(dataIni$DBH,main="DBH",breaks=seq(10,150,20))

Quantiles, diagramme à moustache et fonction de répartition empirique

Remarque : une mesure empirique est une mesure qui décrit l’échantillon, sans référence à la population sous-jacente.

Le quantile empirique d’ordre p d’un échantillon correspond à la plus petite valeur Q présente dans l’échantillon telle qu’au moins une fraction p de l’échantillon a une valeur inférieure ou égale à Q

quantile(dataIni$DBH,0.7,type=1,na.rm=TRUE)  

## 70% 
##  42

#attention, d'autres valeurs de 'type' correspondent à des 
#définitions légèrement différentes de la notion de quantile

Certains quantiles emprique de l’échantillon sont fréquemment utilisés pour sa description : premier décile, premier quartile, médiane, troisième quartile, neuvième décile correspondant respectivement aux quantiles d’ordre 0.1 , 0.25, 0.5, 0.75, 0.9. Les diagrammes à moustache permettent de visualiser la médiane et les quartiles :

boxplot(dataIni$DBH,range=0,ylab="DBH")  

La fonction de répartition est un outil encore plus puissant que le diagramme à moustache, qui permet de visualiser l’ensemble des quantiles de l’échantillon simultanément.

fRep_DBH <- ecdf(dataIni$DBH) 
plot(fRep_DBH,xlab="DBH",ylab="Quantile",main="Fonction de répartition de DBH")

La fonction de répartition de l’échantillon contient toute l’information de l’échantillon. Elle est plus précise que l’histogramme.

Moyenne et variance empiriques

La moyenne empirique de l’échantillon est une valeur située “au milieu” de l’échantillon, elle minimise une certaine distance à l’échantillon. La notion de distance utilisée est la somme des écarts quadratiques avec les valeurs présentes dans l’échantillon. On appellera cette somme coût quadratique par la suite.

coutQuadratique <- function(mu,ech) sum((ech-mu)^2,na.rm=TRUE)
vecteur_m <- 10:140
vecteur_cout <- sapply(vecteur_m,function(m) coutQuadratique(m,dataIni$DBH))
his <- hist(dataIni$DBH,plot=FALSE,breaks=seq(10,140,5))
plot(vecteur_m,vecteur_cout,xlab="Valeur candidate",ylab="Cout quadratique")
polygon(c(his$mids,rev(his$mids)),c(his$counts,rep(0,length(his$counts)))*10000,
        border=NA,col="lightgrey")
points(vecteur_m,vecteur_cout,col=2)
abline(v=mean(dataIni$DBH,na.rm=TRUE),col=4)
text(mean(dataIni$DBH,na.rm=TRUE)+1,max(vecteur_cout),"moyenne",adj=0,col=4)

Remarque : si on avait utilisé une autre notion de distance à l’échantillon, on aurait trouvé une autre quantité que la moyenne de l’échantillon. Par exemple, si on avait minimisé la somme des écarts absolus plutôt que la somme des écarts quadratiques, on aurait obtenu la médiane empirique.

La moyenne des écarts quadratiques obtenus avec la moyenne empirique correspond à la variance empirique de l’échantillon et traduit une notion de dispersion de l’échantillon autour de sa moyenne empirique.

N_DBH <- sum(!is.na(dataIni$DBH))
cout <- coutQuadratique(mean(dataIni$DBH,na.rm=TRUE),dataIni$DBH)/N_DBH
cout #variance empirique de l'échantillon

## [1] 203.7147

Pour remettre dans une unité interprétable la variance empirique, on peut calculer l’écart type empirique.

sqrt(cout) #écart-type empirique de l'échantillon

## [1] 14.27286

Estimer la moyenne et variance de la population

La population d’où provient l’échantillon possède également une moyenne, et une variance, mais celles-ci ne sont pas connues car on n’observe qu’une petite fraction des individus qui la composent, l’échantillon. La moyenne de la population et la moyenne empirique (de l’échantillon) ce n’est donc pas la même chose. Néanmoins, la moyenne empirique permet d’estimer la moyenne de la population, c’est à dire de faire une proposition de valeur dont on peut penser qu’elle n’est pas mauvaise sachant les informations dont on dispose.

Même remarque pour la variance empirique : elle peut permettre d’estimer la variance de la population. Cependant, le fait que la variance empirique soit calculée via des écarts quadratiques de l’échantillon à la moyenne empirique, et non pas via des écarts à la moyenne de la population, mène à une sous-estimation systématique de la variance de la population. Ce problème se traite par l’application d’un facteur correctif à la variance empirique, qui mène à un nouvel estimateur de la variance de la population, qu’on appelle simplement estimateur de la variance dans la suite.

cout*N_DBH/(N_DBH-1) #estimateur de la variance de l'échantillon

## [1] 203.8251

var(dataIni$DBH,na.rm=TRUE) #fonction de R pour calculer directement l'estimateur

## [1] 203.8251

Pour remettre dans une unité interprétable l’estimateur de la variance, on peut considérer sa racine carrée, on parle alors d’estimateur de l’écart type de la population :

sd(dataIni$DBH,na.rm=TRUE)

## [1] 14.27673

Le modèle de population gaussien

Dans la suite, on fait des hypothèses a priori sur la façon dont la variable d’intérêt est distribuée dans la population. On parle de modélisation statistique de la population. Ces hypothèses sont faites pour permettre d’estimer des propriétés de la population tout en quantifiant l’incertitude de l’estimation, et de prédire les futures observations qu’on obtiendra en tirant de nouveaux échantillons dans la population.

Pour illustrer la problématique de l’inférence, rappelons que la moyenne empirique estime la moyenne de la population. On parle d’estimateur de la moyenne. Cependant la moyenne empirique n’est pas exactement égale à la moyenne de la population. On fait une erreur d’estimation, liée au caractère aléatoire de l’échantillonnage. Comment savoir si cette erreur est forte ou faible, comment déterminer à partir de la moyenne empirique les plages de valeur plausible pour la moyenne de la population ? Comment savoir si la moyenne de la population est au dessus ou en dessous d’un certain seuil ? Proposer une modélisation statistique de la population peut permettre d’apporter des réponses à ces questions, sous l’hypothèse que le modèle est valide.

Modélisation gaussienne de la population

La loi normale, loi de Gauss ou loi gaussienne est distribution symétrique complètement définie par la moyenne et la variance, qui ressemble à une cloche. On visualise ci-dessous la densité d’une distribution gaussienne de moyenne \(0\) et variance \(1\) qu’on appelle loi gaussienne centrée réduite.

vX <- seq(-3,3,0.01)
plot(vX,dnorm(vX))

Le modèle gaussien de la population consiste en faire l’hypothèse qu’une variable d’intérêt \(Y\) (par ex. Y=DBH dans notre exemple fil rouge) suit une loi gaussienne au sein de la population et que les échantillons constituent des tirages indépendants dans la population. On note : \[Y_i=μ_Y+E_i\] où \(E_i∼N(0,σ^2_Y)\) désigne les fluctuations des observations autour de la moyenne de la population, \(μ_Y\) est la moyenne de la population et \(σ^2_Y\) est la variance de la population.

Comme dit plus haut la moyenne empirique de l’échantillon est un estimateur de \(μ_Y\), et on la notera \(\hat{\mu}_Y\) dans la suite, la notation “chapeau” signifiant “estimateur de”. De même, on notera \(\hat{\sigma}^2_Y\) l’estimateur de la variance introduit plus haut et \(\hat{\sigma}_Y\) l’estimateur de l’écart-type..

Il existe une fonction pour appliquer le modèle gaussien sur un échantillon dans \(R\):

modGauss <- lm(DBH~1,data=dataIni) #objet qui contient le modèle gaussien de DBH

Évaluer les hypothèses du modèle gaussien

L’écart \(E_i\) entre l’observation \(i\) et la moyenne de la population est l’ erreur asscoiée à l’observation \(i\). Les hypothèses du modèle gaussien portent sur la façon dont les erreurs se distribuent autour de la moyenne de la population : (i) les erreurs sont toutes tirées dans une même loi gaussienne et (ii) les erreurs sont tirées de façon indépendante dans cette loi. L’ensemble des applications qu’on va tirer de ce modèle en termes d’inférence ou de prédiction dépend de la validité de ces hypothèses. Il est donc impératif d’évaluer l’adéquation de ces hypothèses sur la base de l’échantillon dont on dispose avant d’aller plus avant dans la mise en oeuvre. L’évaluation va donc se dérouler en deux étapes : estimer les erreurs associées à chaque élément de l’échantillon, puis sur la base de ces estimateurs, évaluer les hypothèses du modèle gaussien.

Estimer les erreurs du modèle gaussien : notion de résidu

Comme on l’a évoqué plus haut, les erreurs du modèle gaussien sont les écarts entre les valeurs de la variable d’intérêt dans l’échantillon avec la moyenne de la population sous-jacente. Pour l’élément numéro \(i\) de l’échantillon, si on observe la valeur \(y_i\) de la variable d’intérêt, l’erreur associée est \(E_i=y_i-\mu_Y\). Cependant, on ne connaît pas la moyenne de la population \(\mu_Y\) donc on ne peut pas calculer la valeur de cette erreur ! Par contre, on peut estimer la moyenne de la population par la moyenne empirique. Si on utilise la moyenne empirique (\(\hat{\mu}_Y\)) à la place de la moyenne de la population, on obtient une nouvelle quantité : \(\hat{e}_i=y_i-\hat{\mu}_Y\) qu’on peut complètement calculer (c’est l’écart entre l’observation de l’élément \(i\) et la moyenne empirique de l’échantillon). Cet écart \(\hat{e}_i\) est un estimateur de l’erreur du modèle gaussien \(E_i\), on parle de résidu.

Ces résidus peuvent être directement récupérés dans l’objet \(R\) associé au modèle gaussien. Dans l’exemple fil rouge du cours :

modGauss$residuals[1:20] #résidus des 20 premiers éléments de l'échantillon

##         63         64         65         66         67         68         69 
##  12.673713   5.673713   1.673713  23.673713  15.673713   6.673713   7.673713 
##         70         71         72         73         74         75         76 
##  12.673713  15.673713   9.673713   6.673713 -13.326287  25.673713  13.673713 
##         77         78         79         80         81         82 
##  15.673713  19.673713  18.673713   8.673713 -11.326287  11.673713

hist(modGauss$residuals)

Par construction, les résidus ont une moyenne nulle :

mean(modGauss$residuals)

## [1] -5.399894e-15

Dans le cadre du modèle gaussien, quand on divise les résidus empiriques par un estimateur de leur variance, on obtient ce que l’on appelle les résidus standardisés, qui sont censé avoir une variance égale à \(1\).

modGauss_resStand <- rstandard(modGauss)

Adéquation des erreurs à la loi gaussienne

On va évaluer dans quelle mesure les résidus calculés précédemment indiquent que les erreurs sont tirées indépendamment dans une même loi gaussienne.

Si le modèle gaussien est valide et si l’échantillon est grand (typiquement plus de \(30\) individus), les résidus standardisés peuvent être, en première approximation, considérés comme tirés dans une loi gaussienne centrée réduite (de moyenne \(0\) et variance \(1\)). Il s’agit d’une approximation car les résidus empiriques ne sont en pratique ni exactement gaussiens ni indépendants, mais ces écarts deviennent négligeables pour des grands échantillons.

On peut donc évaluer les hypothèses du modèle gaussien en regardant si les quantiles des résidus standardisés sont répartis comme on l’attendrait pour un tirage de \(N\) valeurs dans une loi gaussienne centrée réduite. Pour ce faire, on compare la fonction de répartition des résidus empiriques avec la fonction de répartition théorique d’une loi gaussienne centrée réduite.

Sur l’exemple fil rouge du cours :

modGauss_fRep_ResEmpStand <- ecdf(modGauss_resStand)
modGauss_N <- length(modGauss_resStand)
vP <- pnorm(vX,mean=0,sd=1)
plot(modGauss_fRep_ResEmpStand,xlab="Residus empiriques standardisés",
     ylab="Quantile")
lines(vX,vP,col="blue")
lines(vX,vP+1.96*sqrt(vP*(1-vP)/modGauss_N),lty="dashed",col="blue")
lines(vX,vP-1.96*sqrt(vP*(1-vP)/modGauss_N),lty="dashed",col="blue")

Dans le graphique ci-dessus, la ligne bleu continue montre où on attend les différents quantiles dans une loi normale centrée réduite. En pratique si on tire un échantillon de \(N\) individus dans une loi gaussienne centrée réduite, la fonction de répartition empirique va dévier de cette courbe théorique simplement parce que l’échantillonnage est aléatoire. Ce n’est pas cette déviation “ordinaire” qui nous intéresse, mais la recherche d’une déviation entre la fonction de répartition empirique des résidus et la fonction de répartition théorique de la loi gaussienne centrée réduite qui serait due au fait que les résidus ne viennent pas d’une distribution gaussienne centrée réduite.

On sait prédire mathématiquement la taille-type des fluctuations autour de la fonction de répartition théorique qui sont liées à l’échantillonnage aléatoire. Les fluctuations probables sont encadrées par l’enveloppe des lignes bleues pointillées. La probabilité que la fonction de répartition des résidus empiriques sorte de cette enveloppe à un quantile donné si les hypothèses du modèle gaussien sont correctes est égale à \(0.05\) seulement, autrement dit c’est improbable.

Ici on constate que la fonction de répartition des résidus empiriques standardisés sort de l’enveloppe pour la plupart des quantiles inférieurs à \(0.6\), ce qui suggère que l’hypothèse de distribution gaussienne des résidus du modèle n’est pas en adéquation avec les données.

On peut visualiser différemment cette comparaison des fonctions de répartition avec un diagramme quantile-quantile (QQ-plot en anglais). Ce diagramme consiste en mettre en regard les deux fonctions de répartition, la fonction de répartition théorique de la loi normale, et la fonction de répartition empirique. Pour chaque ordre de quantile, on dessine un point sur le diagramme dont la coordonnée en abscisse (horizontale) vaut la valeur de ce quantile pour la loi normale centrée réduite et la coordonnée en ordonnée (verticale) vaut la valeur de ce quantile pour la fonction de répartition empirique. Par exemple, pour le quantile d’ordre \(0.5\) (la médiane), on dessine un point d’abscisse \(0\) ( médiane de la loi gaussienne centrée réduite) et d’ordonnée égale à la médiane de l’échantillon.

Dans notre exemple fil rouge, le QQ-plot s’obtient comme suit :

library(car)

## Le chargement a nécessité le package : carData

qqPlot(modGauss_resStand,distribution="norm",mean=0,sd=1,line="none")

##  354 1944 
##  115 1541

On a rajouté sur le diagramme l’enveloppe des déviations autour de l’attendu théorique dues à l’échantillonnage aléatoire qu’on attendrait quand les erreurs sont issues d’une loi gaussienne. On retrouve le fait que les quantiles de l’échantillon sortent de l’enveloppe, suggérant que les erreurs ne sont pas distribuées selon la loi gaussienne.

Il est intéressant d’interpréter, à partir du QQ-plot ci-dessus et de l’histogramme des résidus dessiné plus haut, comment la distribution de l’ échantillon dévie de la loi gaussienne centrée réduite. Pour faire simple, on peut distinguer deux types de déviation à la gaussienne : la présence d’une asymétrie (la distribution empirique des échantillons est trop “penchée” d’un côté) et un excès ou défaut de kurtosis (la distribution est trop pointue ou trop carrée au centre).

Lorsque le QQ-plot dévie avec une forme de “banane”, c’est à dire en sortant dans le même sens hors de l’enveloppe aux deux extrêmes de la distribution tout en étant à peu près sur la diagonale pour les valeurs intermédiaires, cela indique un problème d’asymétrie. Une banane avec les pointes vers le haut indique une asymétrie vers la droite, c’est à dire le gros des données légèrement à gauche de zéro et une queue de distribution épaisse vers la droite. C’est le cas pour notre exemple fil rouge. Une banane avec les pointes vers le bas indique une asymétrie dans l’autre sens.

Lorsque le QQ-plot dévie avec une forme en “slalom”, c’est à dire en sortant de l’enveloppe dans des sens différents à chaque extrême tout en étant à peu près sur la diagonale pour les valeurs intermédiaires, cela indique un problème de kurtosis. Si la courbe du QQ-plot est trop faible au début et trop forte à la fin, cela indique un excès de kurtosis, une distribution trop pointue. Si la courbe est au contraire trop forte au début et trop faible à la fin, cela indique un défaut de kurtosis, la distribution empirique des résidus est trop carrée.

On donne ci-dessous quelques exemples de QQ-plot sur des distributions factices pour illustrer l’effet de l’asymétrie et de la kurtosis :

bks <- seq(-8,8,0.1)
vecAlp <- c(1.1,5,1); vecBet <- c(5,1.1,1)
vecTit <- c(
  "Asymétrie à droite",
  "Asymétrie à gauche",
  "Défaut de kurtosis",
  "Excès de kurtosis"
)
par(mfrow=c(length(vecAlp)+1,2))
nul <- sapply(1:length(vecAlp),function(i){
  alp <- vecAlp[i]; bet <- vecBet[i]
  samp <- scale(rbeta(N_DBH,alp,bet))
  hist(samp,breaks=bks,probability=TRUE,main="")
  lines(bks,dnorm(bks,0,1),col="blue")
  title(main=vecTit[i])
  nul <- qqPlot(samp,distribution="norm",mean=0,sd=1,line="none")
})
samp <- scale(rexp(N_DBH,1)*(rbinom(N_DBH,1,0.5)-0.5))
hist(samp,breaks=bks,probability=TRUE,main="")
lines(bks,dnorm(bks,0,1),col="blue")
title(main=vecTit[4])
nul <- qqPlot(samp,distribution="norm",mean=0,sd=1,line="none")

En pratique, les résultats du modèle gaussien (et plus généralement de l’ensemble des modèles linéaires gaussiens présentés dans ce cours) sont robustes aux excès et aux défauts de kurtosis. Ce n’est donc pas très grave d’avoir des QQ-plot en slalom. En revanche l’asymétrie est connue pour générer des biais dans les analyses, il faut donc la traiter.

Transformation Box-Cox

Les transformations de Box-Cox sont une famille de transformations de la variable d’intérêt qui tentent de corriger au mieux les problèmes de déviation des résidus par rapport à une loi gaussienne.

On recherche dans l’exemple fil rouge du cours une transformation de type “puissance” qui permette de régler notre problème d’asymétrie :

modGauss_pT <- powerTransform(modGauss) #recherche d'une transformation Box-Cox de type "puissance"
transfDBH <- dataIni$DBH^modGauss_pT$lambda #application de la tranformation
plot(transfDBH~DBH, data=dataIni,xlab="DBH",ylab="DBH transformé")

On constate que la transformation puissance retenue est décroissante. La variable transformée décroît avec la variable d’origine. On évalue si la transformation corrige le problème d’asymétrie des résidus :

modGauss_BC <- lm(transfDBH~1) 
par(mfrow=c(1,2))
hist(modGauss_BC$residuals/sd(modGauss_BC$residuals),main="",
     xlab="Résidus standardisés")
nul <- qqPlot(modGauss_BC,distribution="norm",line="none")

On a bien une distribution des résidus symétriques (pas de QQ-plot en “banane”) et on a un défaut de kurtosis qui persiste mais qui ne devrait pas trop affecter les inférences et tests présentés ci-dessous.

Indépendance des erreurs

L’hypothèse d’indépendance des erreurs du modèle gaussien revient à postuler qu’ aucun facteur n’a d’influence sur la valeur que prennent les résidus. Pour évaluer ce point, on peut regarder s’il existe un lien entre d’autres covariables du jeu de données et les résidus du modèle gaussien considéré. Pour ce faire, on peut avoir recours à un modèle linéaire gaussien. Nous reviendrons sur ce point dans les sections portant sur le modèle ANOVA et le modèle de régression linéaire.

Dans l’exemple fil rouge du cours, on peut par exemple remarquer que l’échantillonnage des arbres n’est pas aléatoire dans l’espace, mais agrégé en \(9\) relevés :

plot(lati~long, data=dataIni,xlab="Longitude WGS84",ylab="Latitude WGS84",pch=3)
coorPla <- sapply(unique(dataIni$releve),function(idPla){
 vIndPla <- which(dataIni$releve==idPla)
 latiPla <- mean(dataIni$lati[vIndPla],na.rm=TRUE)
 longPla <- mean(dataIni$long[vIndPla],na.rm=TRUE)
 return(c(longPla,latiPla))
})
text(coorPla[1,],coorPla[2,],colnames(coorPla),col=4,cex=0.5,font=2)

Si on visualise la distribution des résidus par relevé à l’aide d’un diagramme à moustache :

vecIndObs <- as.numeric(rownames(modGauss_BC$model))
boxplot(rstandard(modGauss_BC)~dataIni$releve[vecIndObs],xlab="", 
        ylab="Résidu standardisé",las=2,range=0)
abline(h=0,lty="dashed")

Il semblerait que les relevés à gauche du diagramme montrent des résidus avec un biais positif, tandis que les relevés à droite du diagramme (BLO_24 et suivants) montrent des résidus avec un biais négatif. La placette de relevé influencerait donc les résidus, tantôt à la hausse, tantôt à la baisse. En d’autre termes, les résidus tendent à dévier de \(0\) de la même façon dans une même placette, un phénomène qu’on appelle auto-corrélation spatiale qui contredirait l’hypothèse d’indépendance des résidus les uns par rapport aux autres. Néanmoins, il faut considérer ces représentations en diagramme à moustache avec prudence car elles masquent les différences de taille de sous-échantillon entre les différents relevés. On verra une façon plus rigoureuse d’évaluer cet effet d’autocorrélation spatiale dans la section sur le modèle ANOVA. En attendant, on supposera temporairement qu’il n’y a pas de problème de non-indépendance des résidus, par rapport à la placette de relevé ou à toute autre co-variable du jeu de données.

Estimation et test sur la moyenne de la population

On a dit plus haut que la moyenne empirique \(\hat{\mu}_Y\) est un estimateur de la moyenne de la population \(\mu_Y\). Cependant, \(\hat{\mu}_Y\) n’est pas exactement égale à \(\mu_Y\). On fait une erreur d’estimation, liée au caractère aléatoire de l’échantillonnage. L’avantage d’utiliser le modèle gaussien, c’est qu’on connaît la distribution de cette erreur. Plus précisément, on connaît la distribution de la statistique \(t\) :

\[t=\frac{\hat{\mu}_Y−\mu_Y}{\hat{\sigma}_Y/\sqrt{N}}\]

qui suit une loi de Student à \(N-1\) degrés de liberté où \(N\) est la taille de l’échantillon. Cette loi ressemble à une loi gaussienne centrée réduite dès que \(N\) est grand. Grâce à ce résultat on peut produire des intervalles de confiance et faire des tests sur la moyenne la population \(\mu_Y\).

On voit apparaître dans la distribution de l’erreur d’estimation la quantité \(\hat{\sigma}_Y/\sqrt{N}\) qui est nommée erreur standard de la moyenne. Il s’agit de l’écart-type des variations de la moyenne empirique \(\hat{\mu}_Y\) autour de la moyenne de la population \(\mu_Y\).

En résumé, la distribution de \(\hat{\mu}_Y\) autour de \(\mu_Y\) ressemble à une loi gaussienne d’écart-type \(\hat{\sigma}_Y/\sqrt{N}\).

Si on revient à notre exemple fil rouge avec la variable de diamètre transformée, et en supposant les hypothèses du modèle gaussien vérifiées, on peut lire la valeur de la moyenne empirique et de l’erreur standard de la moyenne dans la première et la deuxième colonne de ce tableau :

modGauss_BC_sum <- summary(modGauss_BC)
modGauss_BC_sum$coefficients

##              Estimate   Std. Error  t value Pr(>|t|)
## (Intercept) 0.1628095 0.0007600276 214.2153        0

Le terme “Intercept” désigne ici la moyenne de la population. Les troisième et quatrième colonnes seront interprétées dans la section suivante sur les tests. On peut retrouver manuellement les valeurs de l’estimateur de la moyenne \(\hat{\mu}_Y\) et de l’erreur standard associée :

mean(transfDBH,na.rm=TRUE)

## [1] 0.1628095

sd(transfDBH,na.rm=TRUE)/sqrt(sum(!is.na(transfDBH)))

## [1] 0.0007600276

Intervalle de confiance de la moyenne de la population

Si les hypothèses du modèle gaussien sont vérifiées, la statistique \(t\) présentée ci-dessus suit une loi de Student à \(N-1\) degré de liberté. La valeur de \(t\) a donc par exemple \(95\)% de chances de se trouver entre le quantile \(0.025\) (noté \(q_{t,N-1}(0.025)\)) et le quantile \(0.975\) (noté \(q_{t,N-1}(0.975)\)) de cette loi (on a éliminé les \(2.5\)% les plus extrêmes de la loi de chaque côté). De la même façon, la valeur de \(t\) a \(50\)% de chances de se trouver entre le quantile \(0.25\) (noté \(q_{t,N-1}(0.25)\)) et le quantile \(0.75\) (noté \(q_{t,N-1}(0.75)\)) de cette loi (on a éliminé les \(25\)% les plus extrêmes de la loi de chaque côté).

Si on reste sur l’idée d’un intervalle à \(95\)%, on peut transformer l’encadrement sur \(t\) en un encadrement sur \(\mu_Y\) (que l’on ne connaît pas !) en remarquant que : \[ q_{t,N-1}(0.025) \leq t \leq q_{t,N-1}(0.975) \\ \iff q_{t,N-1}(0.025) \leq \frac{\hat{\mu}_Y−\mu_Y}{\hat{\sigma}_Y/\sqrt{N}} \leq q_{t,N-1}(0.975) \\ \iff q_{t,N-1}(0.025)\hat{\sigma}_Y/\sqrt{N} \leq \hat{\mu}_Y−\mu_Y \leq q_{t,N-1}(0.975) \hat{\sigma}_Y/\sqrt{N} \\ \iff q_{t,N-1}(0.025)\hat{\sigma}_Y/\sqrt{N} -\hat{\mu}_Y\leq −\mu_Y \leq q_{t,N-1}(0.975) \hat{\sigma}_Y/\sqrt{N}-\hat{\mu}_Y \\ \iff \hat{\mu}_Y−q_{t,N-1}(0.025)\hat{\sigma}_Y/\sqrt{N} \geq \mu_Y \geq \hat{\mu}_Y− q_{t,N-1}(0.975) \hat{\sigma}_Y/\sqrt{N} \\ \iff \hat{\mu}_Y−q_{t,N-1}(0.975)\hat{\sigma}_Y/\sqrt{N} \leq \mu_Y \leq \hat{\mu}_Y− q_{t,N-1}(0.025) \hat{\sigma}_Y/\sqrt{N} \]

L’encadrement qu’on obtient sur \(\mu_Y\) est appelé intervalle de confiance, ici à \(95\)%, et il permet de localiser des valeurs plausibles pour la moyenne de la population.

On peut calculer les bornes de l’intervalle de confiance de la moyenne de la population dans le cadre de notre exemple fil rouge avec diamètre transformé :

confint(modGauss_BC,level=0.95)

##                 2.5 %    97.5 %
## (Intercept) 0.1613189 0.1643001

On peut avoir envie de réexprimer ce résultat dans l’échelle naturelle du diamètre des arbres (c’est à dire en cm). Il suffit pour cela d’appliquer sur ces bornes la transformation inverse de la transfromation Box-Cox utilisée précédemment :

confint(modGauss_BC,level=0.95)^(1/modGauss_pT$lambda)

##                2.5 %   97.5 %
## (Intercept) 31.00635 29.95579

On remarque que l’ordre des bornes est inversé, car la transformation Box-Cox retenue ici est décroissante. Cela ne change rien à l’interprétation : le diamètre moyen de la population se trouve avec un niveau de confiance à \(95\)% entre \(30.0\) et \(31.0\)cm (bornes arrondies au mm).

Remarque : on peut avoir aussi envie de réexprimer l’estimateur de la moyenne de la population dans l’échelle naturelle du diamètre des arbres (c’est à dire en cm), en appliquant la même transformation inverse :

mean(transfDBH,na.rm=TRUE)^(1/modGauss_pT$lambda)

##       Y1 
## 30.47414

On a donc une estimation de la moyenne de la population à \(30.5\) cm (arrondie au mm).

Test paramétrique sur la moyenne de la population

Un test paramétrique est une opération visant départager une hypothèse nulle (\(H_0\)) et une hypothèse alternative (\(H_1\)) sur les valeurs d’un ou plusieurs paramètres d’un modèle statistique.

Dans le cas du modèle gaussien, on peut facilement faire des tests sur la valeur de la moyenne de la population \(\mu_Y\) quand l’hypothèse nulle est de la forme “\(\mu_Y\) vaut \(x\)” (par exemple “\(\mu_{DBH}\) vaut \(30\)cm” dans notre cas d’étude). Il y a trois types d’hypothèse alternative (\(H_1\)) classiques : “\(\mu_Y\) est plus grand que \(x\)”, “\(\mu_Y\) est plus petit que \(x\)”, “\(\mu_Y\) est différent de \(x\)”. Dans les deux premiers cas, on parle de test unilatéral, dans le troisième cas on parle de test bilatéral.

On se contentera ici de développer la logique du test bilatéral. L’enjeu du test est de décider si l’on rejette \(H_0\) (“\(\mu_Y\) vaut \(x\)”) au profit de \(H_1\) (“\(\mu_Y\) est différent de \(x\)”) ou non. La décision se prend sur la base de la statistique \(t\) introduite plus haut. Si \(H_0\) est vraie, on peut calculer la valeur de \(t\) puisque : \[t=\frac{\hat{\mu}_Y−x}{\hat{\sigma}_Y/\sqrt{N}}\]

et cette valeur est censée provenir d’un tirage dans une loi de Student à \(N-1\) degrés de liberté. Si \(H_1\) est vraie, la valeur de \(t\) calculée sur la base de l’hypothèse \(H_0\) ne suit plus la loi de Student attendue, car son numérateur \(\hat{\mu}_Y−x\) cumule additivement l’erreur \(\hat{\mu}_Y−\mu_Y\) (qui donne lieu à la loi de Student) et une erreur systématique \(\mu_Y-x\) lié au fait d’avoir fait la mauvaise hypothèse dans le calcul de \(t\). Selon le signe de cette erreur systématique, la statistique \(t\) va tendre à adopter des valeurs trop négatives (si \(\mu_Y<x\)) ou trop positives (si \(\mu_Y>x\)) par rapport à l’attendu de la loi de Student. Dans tous les cas, si \(\mu_Y\neq x\) (\(H_1\)) la statistique \(t\) va donc tendre à prendre des valeurs plus loin de \(0\) qu’attendu sous une loi de Student.

Pour décider si on rejette \(H_0\), on va donc évaluer si \(t\) calculé sous \(H_0\) est “trop” loin de \(0\) par rapport à l’attendu que \(t\) suit une loi de Student à \(N-1\) degrés de liberté. On calcule la p-valeur de ce test, notée \(p\), c’est à dire la probabilité d’observer une valeur de \(t\) encore plus loin de \(0\) que ce qu’on observe alors que \(H_0\) est vraie (et donc que \(t\) suit bien une loi de Student à \(N-1\) degrés de liberté). Si la p-valeur est faible, c’est qu’il est très improbable d’observer une valeur de \(t\) si loin de \(0\) sous \(H_0\), ce qui nous incite à rejeter \(H_0\) au profit de \(H_1\) (“\(\mu_Y\) est différent de \(x\)”). Par convention, on commence à préférer \(H_1\) à \(H_0\) si \(p<0.05\).

L’application du modèle gaussien sous \(R\) contient par défaut un test bilatéral sur la moyenne de la population avec \(H_0 : \mu_Y=0\) et \(H_1 : \mu_Y \neq 0\). En revenant à notre exemple fil rouge, la valeur de la statistique \(t\) sous \(H_0\) et la p-valeur \(p\) du test sont reportées dans les troisième et quatrième colonnes du tableau suivant :

modGauss_BC_sum$coefficients

##              Estimate   Std. Error  t value Pr(>|t|)
## (Intercept) 0.1628095 0.0007600276 214.2153        0

Le résultat du test est sans ambigüité : il est extrêmement improbable d’observer une valeur de \(t\) aussi forte sous \(H_0 (\mu_{DBH}=0)\) ce qui amène à préférer \(H_1 (\mu_{DBH}\neq 0)\), et le fait que \(t>0\) suggère plus précisément que \(\mu_{DBH}> 0\). Cette conclusion est néanmoins triviale puisqu’un diamètre transformé égal à \(0\) correspond à un arbre de diamètre infini. Il est évident que la moyenne des diamètres de la population des arbres n’est pas infinie.

On peut faire des tests plus intéressant fondés sur d’autres valeurs de \(x\). Par exemple, on peut vouloir tester si la moyenne des diamètres des arbres dans la population est significativement différente de \(30\)cm. Cela revient à tester si la moyenne des diamètres transformés est différente de la valeur de \(x\):

x <- 30^modGauss_pT$lambda ; x

##        Y1 
## 0.1641715

ce qu’on peut faire, par exemple en réappliquant le modèle gaussien sur les diamètres transformés translatés de cette valeur et en regardant les colonnes correspondant au test :

modGauss_BC_trans <- lm(I(transfDBH-x)~1)
summary(modGauss_BC_trans)$coefficients

##                 Estimate   Std. Error   t value   Pr(>|t|)
## (Intercept) -0.001361923 0.0007600276 -1.791939 0.07330658

On conclut que la moyenne de la population pour le diamètre transformé est plus faible que \(x\), mais pas de façon significative (\(t\) est négative mais la p-valeur associée est supérieure à \(0.05\)). On en conclut que la moyenne des diamètres dans la population n’est pas significativement différente de \(30\)cm.

Le modèle ANOVA

Motivation : la présence éventuelle de sous-populations

Dans la section précédente sur le modèle gaussien, on a examiné sur notre exemple fil rouge les résidus standardisés en fonction de la placette de relevé :

modGauss_BC_vecIndObs <- as.numeric(rownames(modGauss_BC$model))
boxplot(rstandard(modGauss_BC)~dataIni$releve[modGauss_BC_vecIndObs],xlab="", 
        ylab="Résidu standardisé",las=2,range=0)
abline(h=0,lty="dashed")

Visuellement, on pouvait penser que dans certains relevés, les résidus standardisés des arbres tendaient à être tous au dessus de \(0\) tandis
que dans d’autres les résidus standardisés tendaient à être tous en dessous de \(0\). Cela signifierait que les diamètres transformés des arbres sont tous au dessus de la moyenne de la population dans certains relevés, et tous en dessous de la moyenne de la population dans d’autres. On a souligné que cela contreviendrait à l’hypothèse d’indépendance des erreurs dans le modèle gaussien, on ajoute ici que cela suggère que le population des arbres ne serait en réalité pas homogène, mais plutôt constituée de sous-populations de moyenne différente dans chaque placette.

Le modèle d’analyse de la variance (ANOVA)

On peut modéliser la présence de sous-populations en complexifiant le modèle gaussien de la section précédente. On peut désormais supposer que les individus de l’échantillon sont tirés dans des groupes qu’on sait distinguer (indexés par \(j\) plus bas; des relevés différents dans notre exemple). On permet désormais que ces groupes correspondent à des sous-populations différentes. En d’autres termes, les individus d’un groupe forment un échantillon d’une sous-population pour laquelle la variable à observer \(Y\) est distribuée avec une erreur gaussienne autour d’une moyenne propre à cette sous population.

L’observation faite sur le \(i\)-ème individu du groupe \(j\) est alors modélisée comme suit : \[Y_{ij} = \mu_j +E_{ij}\] où \(\mu_j\) est la moyenne de la sous-population correspondant au groupe \(j\) et \(E_{ij}\) est l’erreur associée à l’observation \(Y_{ij}\). Les erreurs sont toujours supposées indépendantes les unes des autres et distribuées suivant la loi gaussienne \(\mathcal{N}(0,\sigma_Y^2)\). A noter que la variance des erreurs ne change pas entre les sous-populations, on parle d’une hypothèse d’homoscédasticité.

Ce modèle s’apppelle le modèle d’analyse de la variance ou modèle ANOVA. La différence fondamentale avec le modèle gaussien vu plus haut est qu’il y a désormais une valeur moyenne \(\mu_j\) par sous population.

Dans notre exemple, on va appliquer un modèle d’analyse de variance sur le diamètre des arbres en prenant en compte les différents relevés de provenance, afin d’explorer plus avant l’impression visuelle développée en section précédente que les moyennes de diamètre semblent différer entre les relevés.

Le modèle ANOVA s’applique dans R avec la fonction lm déjà utilisée précédemment.

modAnova <- lm(DBH~0+releve,data=dataIni)

Adéquation

Comme dans le modèle gaussien, on commence par vérifier que les hypothèses sur les résidus sont en adéquation avec les données utilisées.

Hypothèse 1 : loi gaussienne des erreurs (ou au moins pas d’asymétrie prononcée)

par(mfrow=c(1,2))
nul <- qqPlot(modAnova$residuals,distribution="norm",line="none")
hist(rstandard(modAnova))

On a un QQ-plot qui suggère un excès de kurtosis mais pas de problème d’asymétrie prononcée. On n’a pas besoin de transformer la variable réponse à ce stade.

Hypothèse 2 : homoscédasticité des erreurs

Pour évaluer l’homoscédasticité des erreurs,on peut utiliser un diagramme de localisation d’échelle (scale-location plot ou SL-plot). Il s’agit d’une représentation de la valeur absolue des résidus empiriques standardisés en fonction de l’estimateur de la moyenne de la sous-population à laquelle ils se rattachent.

Si l’hypothèse d’homoscédasticité est vérifiée, l’amplitude des résidus standardisés doit être la même quelque soit la moyenne de la sous-population à laquelle ils se rattachent, et les points du diagramme doivent être distribués en moyenne autour d’une ligne horizontale. De plus, comme les résidus standardisés sont approximativement distribués comme une loi gaussienne centrée réduite, la ligne horizontale doit se trouver à une valeur de \(0.8\) (ligne bleue dans l’exemple ci-dessous).

On peut évaluer visuellement ce point en effectuant une moyenne lissée sur le SL-plot (ligne rouge dans l’exemple ci-dessous) et en regardant si cette moyenne s’écarte significativement de la ligne horizontale (i.e. est-ce que la ligne bleue sort de l’enveloppe en pointillé).

Dans notre exemple fil rouge, le SL-plot et le diagnostic s’obtiennent comme suit:

#SL-plot
plot(modAnova,which=3,pch=3, add.smooth = FALSE) #Base du SL-plot pré-programmée dans R
abline(h=0.8,col=4,lwd=2) #Ligne horizontale attendue
lo <- loess(sqrt(abs(rstandard(modAnova)))~modAnova$fitted.values) #Moyenne glissante des points
vFit <- sort(unique(modAnova$fitted.values))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~sort(
    unique(modAnova$fitted.values)
  ),col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~sort(
    unique(modAnova$fitted.values)
  ),col=2,lwd=2,lty="dashed")

On constate que les résidus ont une variance trop faible pour les relevés à faible diamètre moyen, et trop forte pour les relevés avec une variance moyenne à forte. Le profil non-monotone de la variance des résidus en fonction des moyennes estimées des sous-populations suggère qu’une transformation monotone de la variable réponse ne résoudra pas le problème.

On cherche d’autres raisons qui expliquent cette fluctuation de la variance. On dispose dans le jeu de donnée d’une information sur les espèces d’arbre correspondant aux observations :

unique(dataIni$recherche_esp_lb_nom_plantae)

##  [1] "Quercus L., 1753"                    "Carpinus betulus L., 1753"          
##  [3] "Fagus sylvatica L., 1753"            "Fraxinus excelsior L., 1753"        
##  [5] "Acer campestre L., 1753"             "Castanea sativa Mill., 1768"        
##  [7] "Prunus avium (L.) L., 1755"          "Acer monspessulanum L., 1753"       
##  [9] "Crataegus monogyna Jacq., 1775"      "Sorbus torminalis (L.) Crantz, 1763"
## [11] "Prunus L., 1753"                     "Sorbus domestica L., 1753"          
## [13] "Quercus canariensis Willd., 1809"    "Castanea Mill., 1754"

Ces arbres ont des tailles très différentes. On peut donc s’attendre à avoir une variance forte des résidus dans des relevés où il y a une grand diversité d’espèces, alors qu’on aura une variance plus faible dans les relevés avec peu d’espèces différentes. Une diversité plus ou moins forte de la composition en espèce entre les différents relevés pourrait donc expliquer l’hétéroscédasticité des résidus.

vRel <- unique(dataIni$releve)
vSpe <- unique(dataIni$recherche_esp_lb_nom_plantae)
tabSpeRel <- sapply(vRel,function(rel){
  iRel <-which(dataIni$releve==rel)
  sapply(vSpe,function(spe){
    sum(dataIni$recherche_esp_lb_nom_plantae[iRel]==spe)
  })
})
vDivSpe <- apply(tabSpeRel,2,function(v) 1-sum(v^2)/sum(v)^2)
plot(sqrt(abs(rstandard(modAnova)))~vDivSpe[match(modAnova$model[,2],vRel)],
     ylab=expression(sqrt(abs(StandardizedResiduals))),xlab="Diversité dans le relevé")
lo <- loess(sqrt(abs(rstandard(modAnova)))~vDivSpe[match(modAnova$model[,2],vRel)])
predLo <- predict(lo,sort(unique(vDivSpe[match(modAnova$model[,2],vRel)])))
lines(predLo~sort(unique(vDivSpe[match(modAnova$model[,2],vRel)])),col=2,lwd=2)

On observe bien une tendance de l’amplitude des résidus à augmenter avec la diversité des espèces dans le relevé, ce qui peut expliquer l’hétéroscédasticité quand on ne prend pas en compte cette diversité d’espèces. Dans la suite de cette section, on va se concentrer sur l’espèce majoritaire du jeu de donnée, le chêne sessile (ici Quercus L., 1753):

dataIni_Querc <- dataIni[
  dataIni$recherche_esp_lb_nom_plantae=="Quercus L., 1753",]
nrow(dataIni_Querc)

## [1] 1467

On réajuste le modèle ANOVA sur ce jeu de donnée réduit, et on refait le diagnostic d’hétéroscédasticité :

modAnova_Querc <- lm(DBH~0+releve,data=dataIni_Querc)
plot(modAnova_Querc,which=3,pch=3, add.smooth = FALSE)
#Ligne horizontale attendue
abline(h=0.8,col=4,lwd=2)
#Moyenne glissante
lo <- loess(sqrt(abs(rstandard(modAnova_Querc)))~modAnova_Querc$fitted.values)
vFit <- sort(unique(modAnova_Querc$fitted.values))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")

On a toujours un problème d’hétroscédasticité. Cependant, contrairement au cas précédent, la variance des résidus augmente de façon monotone avec l’estimateur de la moyenne de la sous population correspondante, ce qui suggère qu’une transformation concave peut régler le problème. On tente donc une transformation Box-Cox:

modAnova_Querc_pT <- powerTransform(modAnova_Querc)
modAnova_Querc_pT$lambda

##         Y1 
## -0.1123729

modAnova_Querc_pT$roundlam

##         Y1 
## -0.1123729

L’exposant de transformation suggéré est proche de 0. Dans ce cas, on peut aussi bien appliquer une transformation de type logarithme, plus simple à interpréter :

modAnova_Querc_BC <- lm(I(log10(DBH))~0+releve,data=dataIni_Querc)

On examine pour la troisième fois l’homoscédasticité :

plot(modAnova_Querc_BC,which=3,pch=3, add.smooth = FALSE)
#Ligne horizontale attendue
abline(h=0.8,col=4,lwd=2)
#Moyenne glissante
lo <- loess(
  sqrt(abs(rstandard(modAnova_Querc_BC)))~modAnova_Querc_BC$fitted.values) 
vFit <- sort(unique(modAnova_Querc_BC$fitted.values))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")

Le résultat est beaucoup mieux ! On vérifie qu’on a toujours une distribution symétrique et approximativement gaussienne des erreurs :

qqPlot(modAnova_Querc_BC,distribution="norm",line ="none")

## 370 987 
## 117 419

On a préservé une distribution des résidus en adéquation avec l’hypothèse de distribution gaussienne des erreurs (moyennant un léger excès de kurtosis sans conséquence importante).

Hypothèse 3 : indépendance des erreurs

Maintenant qu’on a filtré pour ne préserver qu’une espèce et qu’on intègre l’effet du relevé dans le modèle, on n’a plus de variable disponible pour laquelle on peut penser qu’elle génèrerait une dépendance dans les résidus, sauf peut-être la position précise des arbres au sein des relevés, mais cela dépasse le cadre de ce cours.

Estimation et tests sur les moyennes des sous-populations

Estimation et intervalle de confiance des moyennes des sous-populations

Dans le modèle ANOVA, l’estimateur \(\hat{\mu}_j\) de la moyenne de la sous-population \(j\) correspond à la moyenne empirique des échantillons appartenant au groupe \(j\). En utilisant le même type d’arguments que dans le modèle gaussien, mais appliqué cette fois aux sous-populations, on peut obtenir des valeurs d’erreur strandard de la moyenne associées à ces sous-populations. Ces informations sont résumées dans les deux premières colonnes du tableau suivant :

summary(modAnova_Querc_BC)$coefficients

##              Estimate  Std. Error  t value Pr(>|t|)
## releveBLO_1  1.657372 0.007354809 225.3453        0
## releveBLO_12 1.371677 0.004579467 299.5275        0
## releveBLO_13 1.779497 0.010570408 168.3470        0
## releveBLO_17 1.394012 0.004425548 314.9919        0
## releveBLO_21 1.661432 0.007158645 232.0874        0
## releveBLO_24 1.673134 0.009597117 174.3371        0
## releveBLO_27 1.671012 0.007836013 213.2477        0
## releveBLO_4  1.603949 0.007665274 209.2488        0
## releveBLO_9  1.653019 0.006369451 259.5229        0

On peut remettre les estimateurs des moyennes des sous-populations dans leur échelle naturelle en appliquant la transformation inverse de celle ayant servi à transformer la variable réponse :

10^(summary(modAnova_Querc_BC)$coefficients[,1])

##  releveBLO_1 releveBLO_12 releveBLO_13 releveBLO_17 releveBLO_21 releveBLO_24 
##     45.43304     23.53296     60.18619     24.77488     45.85973     47.11227 
## releveBLO_27  releveBLO_4  releveBLO_9 
##     46.88264     40.17439     44.97992

On peut également générer des intervalles de confiance sur la moyenne des sous- populations :

confint(modAnova_Querc_BC,level=0.95)

##                 2.5 %   97.5 %
## releveBLO_1  1.642945 1.671799
## releveBLO_12 1.362694 1.380660
## releveBLO_13 1.758762 1.800232
## releveBLO_17 1.385330 1.402693
## releveBLO_21 1.647389 1.675474
## releveBLO_24 1.654308 1.691960
## releveBLO_27 1.655641 1.686383
## releveBLO_4  1.588913 1.618985
## releveBLO_9  1.640524 1.665513

Ici encore, on peut remettre ces intervalles de confiance sur l’échelle naturelle du paramètre :

tabConfInt <- 10^confint(modAnova_Querc_BC,level=0.95); tabConfInt

##                 2.5 %   97.5 %
## releveBLO_1  43.94856 46.96767
## releveBLO_12 23.05120 24.02480
## releveBLO_13 57.38019 63.12941
## releveBLO_17 24.28457 25.27509
## releveBLO_21 44.40063 47.36678
## releveBLO_24 45.11369 49.19939
## releveBLO_27 45.25233 48.57169
## releveBLO_4  38.80727 41.58967
## releveBLO_9  43.70432 46.29275

On peut restituer cette analyse sur une figure :

x <- barplot(10^(summary(modAnova_Querc_BC)$coefficients[,1]),las=2)
arrows(x0 = x, y0= tabConfInt[,1],
  y1=tabConfInt[,2], code=3,angle=90,
  length=0.05)

Tests

Comparer les moyennes des sous-populations à une valeur fixe

Les tests du modèle ANOVA sur les moyennes des sous-populations se déroulent comme le test sur la moyenne de la population dans le modèle gaussien. Par défaut, les troisième et quatrième colonnes du tableau ci-dessous reportent les résultats (statistique de test \(t\) et p-valeur) associés à un test bilatéral de l’hypothèse nulle \(H_0 : \mu_j=0\) contre l’hypothèse alternative \(H_1 : \mu_j \neq 0\).

summary(modAnova_Querc_BC)$coefficients

##              Estimate  Std. Error  t value Pr(>|t|)
## releveBLO_1  1.657372 0.007354809 225.3453        0
## releveBLO_12 1.371677 0.004579467 299.5275        0
## releveBLO_13 1.779497 0.010570408 168.3470        0
## releveBLO_17 1.394012 0.004425548 314.9919        0
## releveBLO_21 1.661432 0.007158645 232.0874        0
## releveBLO_24 1.673134 0.009597117 174.3371        0
## releveBLO_27 1.671012 0.007836013 213.2477        0
## releveBLO_4  1.603949 0.007665274 209.2488        0
## releveBLO_9  1.653019 0.006369451 259.5229        0

Ici l’hypothèse \(\mu_j=0\) correspond dans l’échelle naturelle à une moyenne de diamètre dans la sous-population \(j\) de l’ordre de \(1\)cm, ce qui est invraisemblalement faible. Il n’est donc pas surprenant de retreouver des valeurs de \(t\) positives avec des p-valeurs très faibles dans toutes les sous-populations : les moyennes sont très probablement supérieures à \(1\)cm dans toutes les sous-populations.

On peut faire un test sur les moyennes des sous populations par rapport à une valeur plus intéressante, par exemple \(25\)cm. Dans l’échelle transformée, cela correpond à une valeur de de \(\log_{10}(25)=1.4\). Pour faire ce test, on réapplique le modèle ANOVA sur une variable réponse translatée de la quantité par rapport à laquelle on veut tester:

summary(lm(I(log10(DBH)-log10(25))~0+releve,data=dataIni_Querc))$coefficients

##                  Estimate  Std. Error    t value      Pr(>|t|)
## releveBLO_1   0.259431817 0.007354809 35.2737680 2.048468e-197
## releveBLO_12 -0.026263369 0.004579467 -5.7350269  1.184205e-08
## releveBLO_13  0.381556860 0.010570408 36.0967023 3.231512e-204
## releveBLO_17 -0.003928474 0.004425548 -0.8876809  3.748594e-01
## releveBLO_21  0.263491494 0.007158645 36.8074528 4.245112e-210
## releveBLO_24  0.275194028 0.009597117 28.6746555 1.115559e-143
## releveBLO_27  0.273072084 0.007836013 34.8483434 6.657882e-194
## releveBLO_4   0.206009295 0.007665274 26.8756610 1.728487e-129
## releveBLO_9   0.255078636 0.006369451 40.0471920 6.343245e-237

Pour certains relevés la moyenne du diamètre transformé pour la sous-population est significativement plus grande que \(\log_{10}(25)\) (par exemple BLO_1), pour certains relevés la moyenne du diamètre transformé pour la sous-population n’est pas différente de \(\log_{10}(25)\) de façon significative (par exemple BLO_17 où la p-valeur est de l’ordre de \(0.37\)), pour certains relevés la moyenne du diamètre transformé pour la sous-population est significativement plus petite que \(\log_{10}(25)\) (par exemple BLO_12).

Comparer les moyennes des sous populations à une sous-population de référence

On peut souhaiter comparer les moyennes des sous-populations entre elles, plutôt que de les comparer à une valeur fixe. Par exemple, on peut souhaiter savoir si les sous-populations ont une moyenne plus ou moins forte qu’une sous-population de référence qu’on désignera comme la population \(1\) par convention. Dans ce cas, il faut reformuler légèrement le modèle ANOVA :

\[Y_{ij} = \mu_j +E_{ij} \\ \iff Y_{ij} = \mu_1+\mu_j-\mu_1 +E_{ij} \\ \iff Y_{ij} = \mu_1+\delta_{1j} +E_{ij} \]

où \(\mu_1\) est la moyenne de la sous-population \(1\) et \(\delta_{1j}=\mu_j-\mu_1\) est l’écart entre la moyenne de la sous-population \(j\) et la moyenne de la sous-population \(1\).

On insiste qu’il s’agit exactement du même modèle ANOVA que précédemment, mais avec un jeu de réécriture pour faire apparaître les paramètres \(\delta_{1j}\) pour lesquels il est intéressant de tester si leur valeur est différente de \(0\).

Pour appliquer cette réécriture dans \(R\) il suffit de modifier légèrement la définition du modèle. Dans notre exemple :

#Définition d'une sous-population de référence, ici BLO_11 :
dataIni_Querc$releve <- relevel(as.factor(dataIni_Querc$releve),ref="BLO_1")
#Reformulation du modèle ANOVA :
modAnovaDif_Querc_BC <- lm(I(log10(DBH))~releve,data=dataIni_Querc)
summary(modAnovaDif_Querc_BC)$coefficients

##                  Estimate  Std. Error     t value      Pr(>|t|)
## (Intercept)   1.657371826 0.007354809 225.3453329  0.000000e+00
## releveBLO_12 -0.285695186 0.008663991 -32.9750109 1.718816e-178
## releveBLO_13  0.122125043 0.012877373   9.4836924  9.628104e-21
## releveBLO_17 -0.263360291 0.008583629 -30.6816959 8.452238e-160
## releveBLO_21  0.004059677 0.010263499   0.3955451  6.924987e-01
## releveBLO_24  0.015762211 0.012091231   1.3036068  1.925742e-01
## releveBLO_27  0.013640267 0.010746921   1.2692255  2.045639e-01
## releveBLO_4  -0.053422522 0.010623071  -5.0289152  5.546425e-07
## releveBLO_9  -0.004353182 0.009729497  -0.4474210  6.546378e-01

Le paramètre “Intercept” désigne la moyenne de la sous-population de référence. Les autres paramètres sont maintenant les \(\delta_{1j}\), les écarts entre la moyenne des sous-population et cette référence. On a comme précédemment un estimateur (première colonne), une erreur standard (deuxième colonne), une statistique \(t\) sous l’hypothèse nulle \(\delta_{1j}=0\) (troisième colonne) et la p-valeur du test bilatéral associé (quatrième colonne). On peut voir que la sous-population du relevé BLO_21 est a une moyenne de diamètre transformé significativement plus forte que la sous population de référence BLO_11 (statistique \(t\) positive et p_valeur environ égale à \(0.01\)).

Coefficients de détermination

Si l’adéquation des hypothèses du modèle ANOVA aux données est satisfaisante, on peut évaluer ce qu’on a gagné à poser l’hypothèse qu’il y a des sous-populations distinctes associées à chaque groupe. C’est le rôle que jouent les coefficients de détermination.

Coefficient de détermination \(R^2\)

Dans le modèle gaussien, on a vu qu’on obtient un estimateur de la moyenne de la population en considérant la moyenne empirique, c’est à dire le nombre qui minimise la distance à l’échantillon qu’on a appelé coût quadratique.

Le modèle ANOVA postule l’existence d’une sous population sous-jacente à chaque groupe de l’échantillon. Dans chaque sous population, la moyenne est estimée en prenant la moyenne empirique du sous-échantillon correspondant. Si on appelle \(J\) le nombre de groupes, on a donc \(J\) estimateurs de moyennes, un par sous-population. Avec ces \(J\) estimateurs, on définit un nouveau coût quadratique, égal à la somme des coûts de chacune de ces moyennes empiriques dans son sous-échantillon.

Ce nouveau coût est plus faible que celui qu’on avait avec une seule moyenne globale, puisque chaque moyenne empirique est plus près de son sous échantillon que ne l’était la moyenne globale. Dans notre exemple, on peut visualiser cette idée avec le graphique suivant :

x <- boxplot(log10(DBH)~releve,data=dataIni_Querc,range=0,las=2,xlab="",cex.axis=0.5)
vCoef <- modAnova_Querc_BC$coefficients[
    match(x$names,gsub("releve","",names(modAnova_Querc_BC$coefficients)))]
points(1:length(vCoef),rep(mean(log10(dataIni_Querc$DBH),na.rm=TRUE),length(vCoef)),pch=21,bg=4)
points(1:length(vCoef),
  modAnova_Querc_BC$coefficients[
    match(x$names,gsub("releve","",names(modAnova_Querc_BC$coefficients)))],
  pch=21,bg=2)

Pour chaque diagramme à moustache, qui correspond à un sous échantillon pour un relevé, le point rouge, qui est la moyenne empirique pour le sous-échantillon considéré, est plus “près” du sous échantillon que le point bleu, qui est la moyenne empirique globale, sur tout l’échantillon. Mathématiquement, le coût quadratique des points rouges (somme des coûts quadratiques des points rouges dans chaque sous échantillon) est plus faible que le coût quadratique associé à la moyenne globale (somme des coûts quadratiques des points bleus) :

modGauss_Querc_BC <- lm(I(log10(DBH))~1,data=dataIni_Querc)
barplot(c(
    sum(modGauss_Querc_BC$residuals^2),
    sum(modAnova_Querc_BC$residuals^2)),
  col=c(4,2),
  names.arg=c("sans relevé","avec relevé"),ylab="Coût quadratique")

On quantifie cette diminution du coût quadratique induite par le fait de distinguer les sous-échantillons associés à chaque relevé par un indice entre \(0\) et \(1\) noté \(R^2\), et appelé coefficient de détermination. Ce dernier se calcule comme suit :

\(R^2 = 1- \frac{\text{coût quadratique avec sous-populations}}{ \text{coût quadratique sans sous-populations}}\)

Le \(R^2\) quantifie donc la fraction du coût quadratique qui a été absorbée par l’utilisation d’une moyenne par sous population plutôt qu’une moyenne globale.

Le \(R^2\) figure parmi les sorties classiques de la fonction lm du logiciel R quand le modèle ANOVA est écrit en termes de différences entre sous-populations :

summary(modAnovaDif_Querc_BC)$r.squared

## [1] 0.7493082

Avec le modèle ANOVA écrit directement en termes de moyennes des sous-populations, il faut le calculer à la main :

1-sum(modAnova_Querc_BC$residuals^2)/sum(modGauss_Querc_BC$residuals^2)

## [1] 0.7493082

Pour les analyses de coefficient de détermination il est donc plus simple (et recommandé) d’écrire le modèle ANOVA sous forme de différence entre les sous-populations et une population de référence.

Coefficient de détermination \(R^2\) ajusté

Le coeffcient de détermination ajusté quantifie de combien on a réduit notre incertitude sur la valeur que prendrait un nouvel individu tiré dans une sous-population avec le modèle ANOVA par rapport au modèle gaussien sans sous-populations.

\(R^2_{adj} = 1- \frac{\text{estimateur de la variance de l'erreur dans les sous-populations}}{\text{estimateur de la variance de l'erreur dans la population totale}}\)

summary(modAnovaDif_Querc_BC)$adj.r.squared

## [1] 0.7479279

Tandis que le \(R^2\) se focalise sur l’amélioration de la qualité de description de l’échantillon (via la baisse du coût quadratique), le \(R^2_{adj}\) se focalise donc sur l’amélioration des capacités prédictives. Le \(R^2_{adj}\) est toujours plus faible que le \(R^2\). Contrairement au \(R^2\) qui est mathématiquement compris entre \(0\) et \(1\), le \(R^2_{adj}\) peut dans certains cas devenir négatif : quand la variance dans les sous-populations est mal estimée car les sous-échantillons associés sont petits, alors que la variance globale de la population reste bien estimée, car il y a beaucoup de groupes, ou dit autrement, quand il y a trop de groupes par rapport à la taille de l’échantillon.

Test de Fisher ou \(F\)-test

N’importe quel découpage d’un échantillon en sous-groupe mène à un modèle ANOVA avec un \(R^2>0\), même si ces sous-groupes sont arbitraires et ne correspondent pas à des sous-populations avec des moyennes différentes. Il est donc important de disposer d’un moyen d’évaluer si le \(R^2\) obtenu est supérieur à ce qu’on obtiendrait par pur hasard, simplement su fait de l’échantillonnage aléatoire dans chaque sous-échantillon. On dispose pour cela d’un test statistique, le F-test, qui considère les hypothèses nulle et aletrnative suivantes :

• \(H_0 : \forall j, \mu_{j}=\mu_1\) (toutes les sous-populations ont la même moyenne) ;
• \(H_1 : \exists j, \mu_j\neq\mu_1\) (au moins une des sous-populations a une moyenne qui diffère des autres)

L’hypothèse \(H_0\) revient à dire que le modèle ANOVA ne sert à rien et qu’un modèle gaussien aurait suffi. Le F-test se pratique sur une statistique \(F\) dérivée du \(R^2\) qui se calcule comme suit :

\[F=\frac{\text{df}_\text{res}}{\text{df}_\text{mod}-1}\frac{R^2}{1-R^2}\]

Dans le calcul de \(F\), \(\text{df}_\text{mod}\) correspond au nombre de degrés de liberté consommés par le modèle ANOVA par rapport au modèle gaussien. Il correspond du nombre de sous-populations introduites dans le modèle ANOVA. Pour quoi appelle-t-on “consommation d’un degré de liberté” le fait d’introduire une sous-population ? Parce que qu’au sein de l’échantillon associé à cette sous-population, on estime la moyenne de la population par la moyenne empirique, et cette moyenne devient une connaissance du modèle. Sachant cette moyenne, les éléments du sous-échantillon ne sont plus libres de varier indépendamment les uns des autres : si on les connaît tous sauf un, on peut déduire le dernier grâce à la moyenne. On dit qu’on a perdu un degré de liberté dans le sous-échantillon, en le faisant basculer dans la connaissance du modèle.

Avant application du modèle, il y a un degré de liberté par élément de l’échantillon. La quantité \(\text{df}_\text{res}\) correspond au nombre de degrés de liberté qu’il reste dans l’échantillon une fois déduits ceux consommés par le modèle. Mathématiquement parlant on a : \(\text{df}_\text{res}=N-\text{df}_\text{mod}\).

Sous l’hypothèse \(H_0\), on connaît la distribution de la statistique \(F\) :

\[F \sim \mathcal{F}(\text{df}_\text{mod}-1,\text{df}_\text{res})\]

où \(\mathcal{F}\) est une loi de Fisher (qui donne son nom au test). Cette loi a deux paramètres, qui dépendent des degrés de liberté consommés par le modèle et résiduels.

On va rejeter \(H_0\) si \(F\) prend une valeur trop forte par rapport à l’attendu donné par la loi de Fisher. En effet une valeur forte de \(F\) signifie que le coefficient de détermination \(R^2\) est plus fort qu’attendu sous \(H_0\), donc que le modèle ANOVA décrit mieux l’échantillon que ce que l’on aurait attendu si les moyennes des sous-populations étaient toutes identiques et qu’un modèle gaussien avait suffi.

Dans notre exemple, on obtient le informations relatives à la statistique \(F\) comme suit:

summary(modAnovaDif_Querc_BC)$fstatistic

##     value     numdf     dendf 
##  542.8701    8.0000 1453.0000

où on peut lire dans l’ordre : la valeur de la statistique \(F\), la valeur de \(\text{df}_\text{mod}-1\) et la valeur de \(\text{df}_\text{res}\). Attention au fait que les calculs sont effectués sur la formulation différentielle du modèle ANOVA (le modèle avec les \(\delta_{1j}\)).

La p-valeur du test (rejet ou non de \(H_0\)) fondé sur la statistique \(F\) se lit en fin de résumé du modèle :

summary(modAnovaDif_Querc_BC)

## 
## Call:
## lm(formula = I(log10(DBH)) ~ releve, data = dataIni_Querc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41777 -0.04946  0.00019  0.05315  0.32729 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   1.657372   0.007355 225.345  < 2e-16 ***
## releveBLO_12 -0.285695   0.008664 -32.975  < 2e-16 ***
## releveBLO_13  0.122125   0.012877   9.484  < 2e-16 ***
## releveBLO_17 -0.263360   0.008584 -30.682  < 2e-16 ***
## releveBLO_21  0.004060   0.010263   0.396    0.692    
## releveBLO_24  0.015762   0.012091   1.304    0.193    
## releveBLO_27  0.013640   0.010747   1.269    0.205    
## releveBLO_4  -0.053423   0.010623  -5.029 5.55e-07 ***
## releveBLO_9  -0.004353   0.009729  -0.447    0.655    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08256 on 1453 degrees of freedom
##   (5 observations effacées parce que manquantes)
## Multiple R-squared:  0.7493, Adjusted R-squared:  0.7479 
## F-statistic: 542.9 on 8 and 1453 DF,  p-value: < 2.2e-16

Dans l’exempple, la p-valeur obtenue est très proche de \(0\), ce qui signifie qu’il est peu vraisemblavle d’observer une valeur de \(F\) (et donc un \(R^2\)) aussi fort s’il n’y a pas de différence de moyenne entre les sous-populations. On rejette donc \(H_0\), et on conclut qu’il était pertinent d’introduire des groupes fondés sur les relevés pour analyser le jeu de données.

Remarque : l’analyse du modèle ANOVA fondé sur les relevés menée dans cette section est une façon rigoureuse de tester l’impression visuelle qu’on avait dans la section sur le modèle modèle gaussien sur la non-indépendance des résidus.

Modèle de régression linéaire

Dans le modèle ANOVA, on dispose de sous échantillons bien délimités, identifiés par des étiquettes. Cependant, on peut avoir besoin de définir des sous populations sur la base d’un critère plus continu. Par exemple, on supposer que l’effet du relevé qu’on avait plus haut est en fait un effet de l’altitude, qui jouerait sur le diamètre des arbres. Les relevés se situant à des altitudes différentes, les moyennes des populations d’arbres sont différentes entre les relevés, mais des arbres à des altitudes différentes au sein d’un même relevé en pente pourraient également correspondre à des populations différentes.

Si on suit cette logique au bout, on se une sous-population par valeur du critère continu : la population des arbres à une altitude de \(100\)m, la population des arbres à une altitude de \(101\)m etc (voire des catégories plus fines si la mesure de l’altitude le permet). Il y a alors très peu d’observations (souvent une seule) par sous-population. Cette situation est inexploitable telle quelle par une approche de type ANOVA : on a déjà évoqué que le \(R^2_{adj}\) sera négatif, le \(R^2\) vaudra \(1\), aucun test ne sera possible et l’intérêt du modèle sera nul, car il sera aussi compliqué que le jeu de donnée initial.

Une option est bien sûr de créer des catégories fondées sur des plages du critère continu, et de revenir sur un modèle ANOVA qui compare ces différentes catégories, mais il faut alors de l’expertise pour créer un bon découpage selon l’objet d’étude. Si on fait plus d’hypothèses sur la façon dont la variable d’intérêt dépend du critère, d’autres stratégies sont possibles dont la régression linéaire.

Dans le cadre de la régression linéaire, on suppose que les moyennes des sous populations ne varient pas arbitrairement avec la valeur du critère continu, mais de façon linéaire. Formellement, le modèle devient: \[Y_{i} = \mu + \beta x_i + E_{i}\] où \(Y_{i}\) est la valeur de la variable d’intérêt pour l’observation \(i\), \(x_i\) est la valeur pour l’observation \(i\) de la covariable qui sert de critère. Le paramètre \(\mu\) est l’ordonnée à l’origine (\(x=0\)) de la relation linéaire, et \(\beta\) est le coeffciient directeur de la relation. L’erreur \(E_{i}\) décrit la fluctuation des invidus de la population \(x=x_i\) autour de la moyenne \(\mu+\beta x_i\). On conserve les trois mêmes hypothèses sur les erreurs que pour le modèle ANOVA :

• les erreurs \(E_{i}\) sont distribuées suivant une loi gaussienne ;
• les erreurs \(E_{i}\) sont homoscédastiques (même variance \(\sigma^2\)) ;
• les erreurs \(E_{i}\) sont indépendantes les unes des autres.

L’application du modèle linéaire se fait en utilisant la fonction lm de R, comme pour le modèle ANOVA :

modReg_Querc <- lm(DBH~alti,data=dataIni_Querc)

Adéquation des hypothèses aux données

Distribution gaussienne des erreurs

qqPlot(modReg_Querc,distribution="norm",line="none")

##  354 1986 
##  104 1311

On trouve un QQ-plot en banane symptomatique d’une distribution asymétrique des résidus. On effectue une transformation des données Box-Cox pour tenter d’atténuer le problème.

modReg_Querc_pT <- powerTransform(modReg_Querc) #roundLam =0 --> on utilise le log
modReg_Querc_BC <- lm(I(log10(DBH))~alti,data=dataIni_Querc)
qqPlot(modReg_Querc_BC,distribution="norm",line="none")

##  354 1986 
##  104 1311

Le diagramme du QQ-plot ne témoigne plus d’une asymétrie prononcée.

Homoscédasticité des erreurs

plot(modReg_Querc_BC,which=3,pch=3, add.smooth = FALSE)
#Ligne horizontale attendue
abline(h=0.8,col=4,lwd=2)
#Moyenne glissante
lo <- loess(
  sqrt(abs(rstandard(modReg_Querc_BC)))~modReg_Querc_BC$fitted.values) 
vFit <- sort(unique(modReg_Querc_BC$fitted.values))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")

On constate une légère hétéroscédasticité sur le SL-plot : les sous-populations de moyennes intermédiaire à forte tendent à avoir une variance des résidus standardisés supérieure à celle des sous populations de moyenne faible.

Indépendance des erreurs

On fait un focus sur un troisième diagramme de diagnostic qui permet de visualiser simplement la valeur des résidus standardisés en fonction de la moyenne de la sous population :

plot(modReg_Querc_BC$residuals~modReg_Querc_BC$fitted.values)
abline(h=0,col=4)
#Moyenne glissante
lo <- loess(
 modReg_Querc_BC$residuals~modReg_Querc_BC$fitted.values) 
vFit <- sort(unique(modReg_Querc_BC$fitted.values))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")

On constate que les résidus ont tendance à montrer des anomalies négatives au centre et positives sur les bords. S’ils se comportaient de façon indépendante les uns des autres, on n’aurait pas cette structure.

On peut regarder le même diagramme en fonction de l’altitude.

plot(rstandard(modReg_Querc_BC)~modReg_Querc_BC$model$alti)
prof <- loess(rstandard(modReg_Querc_BC)~modReg_Querc_BC$model$alti)
yProf <- predict(prof); xProf <- modReg_Querc_BC$model$alti
lines(yProf[order(xProf)]~sort(xProf),col=2)
abline(h=0,col="grey")

On constate donc un problème de non indépendance des résidus en fonctions de l’altitude, ce qui suggère que la relation linéaire qu’on a postulé à l’origine ne capte pas bien le profil de variation du log-diamètre en fonction de l’altitude. Cette relation doit être courbe. On peut visualiser la moyenne de log-diamètre prédite (en bleu ci-dessous) par le modèle en fonction de l’altitude, et superposer les log-diamètres observés pour se rendre compte du phénomène :

plot(modReg_Querc_BC$model[,1]~modReg_Querc_BC$model[,2],xlab="Altitude",ylab=expression(log[10](DBH)))
lines(sort(modReg_Querc_BC$model$alti),modReg_Querc_BC$coefficients[1] + modReg_Querc_BC$coefficients[2]*sort(modReg_Querc_BC$model$alti),col=4)

Le comportement non indépendant des résidus en fonction de l’altitude suggère d’introduire un effet quadratique de l’altitude. Il s’agit d’un modèle du type

\[Y_{i} = \mu + \beta x_i + \gamma x_i^2 + E_{i}\]

qui permet de capter une relation courbe. Il s’applique comme suit :

modRegQuad_Querc <- lm(DBH~alti+I(alti^2),data=dataIni_Querc)

et s’analyse de la même façon que les modèles précédents. On déroule les diagrammes de diganostic ci-dessous

#DISTRIBUTION SYMETRIQUE
qqPlot(modRegQuad_Querc,distribution="norm",line="none")

##  354 1986 
##  104 1311

modRegQuad_Querc_pT <- powerTransform(modRegQuad_Querc)
modRegQuad_Querc_BC <- lm(I(log10(DBH))~alti+I(alti^2),data=dataIni_Querc)
qqPlot(modRegQuad_Querc_BC,distribution="norm",line="none")

##  987 1899 
##  419 1242

#HOMOSCEDASTICITE
plot(modRegQuad_Querc_BC,which=3,pch=3, add.smooth = FALSE)
#Ligne horizontale attendue
abline(h=0.8,col=4,lwd=2)
#Moyenne glissante
lo <- loess(
  sqrt(abs(rstandard(modRegQuad_Querc_BC)))~modRegQuad_Querc_BC$fitted.values) 
vFit <- sort(unique(modRegQuad_Querc_BC$fitted.values))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")

#INDEPENDANCE
plot(modRegQuad_Querc_BC$residuals~modRegQuad_Querc_BC$fitted.values)
abline(h=0,col=4)
#Moyenne glissante
lo <- loess(
 modRegQuad_Querc_BC$residuals~modRegQuad_Querc_BC$fitted.values) 
vFit <- sort(unique(modRegQuad_Querc_BC$fitted.values))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")

On constate qu’on a atténué les problèmes d’hétroscédasticité et de non indépendance en fonction de l’altitude sans les éliminer complètement, on considèrera donc les résultats ci-dessous avec prudence.

Estimation, coefficient de détermination et F-test

summary(modRegQuad_Querc_BC)

## 
## Call:
## lm(formula = I(log10(DBH)) ~ alti + I(alti^2), data = dataIni_Querc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.3897 -0.1077 -0.0154  0.1131  0.3907 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.640e+00  8.742e-02   41.63   <2e-16 ***
## alti        -1.235e-02  4.936e-04  -25.01   <2e-16 ***
## I(alti^2)    1.744e-05  6.775e-07   25.74   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1357 on 1458 degrees of freedom
##   (6 observations effacées parce que manquantes)
## Multiple R-squared:  0.3199, Adjusted R-squared:  0.319 
## F-statistic:   343 on 2 and 1458 DF,  p-value: < 2.2e-16

Les éléments du résumé s’interprète de façon similaire au modèle ANOVA. Concernant les estimateurs des termes du modèles l’estimation de l’ordonnée à l’origine (Intercept) vaut \(2.810e+00\), l’estiumation du coefficient \(\beta\) (alti) vaut \(-6.532e-03\) et l’estimation du coefficient gamma (I(alti^2)) vaut \(9.299e-06\). On a chaque fois une erreur standard qui permet d’apprécier l’intervalle de confiance associé à cet estimateur, et on a pour chaque paramètre un test de Student qui évalue s’il est différent de 0. Ici on juge que les trois paramètres sont significativement différents de \(0\). On a une information sur le fait que \(6\) observations ont été supprimée car l’altitude n’était pas renseignée. On a ensuite les coefficients de détermination \(R^2\) (Multiple R-squared) et \(R^2_{adj}\) (Adjusted R-squared) puis le test \(F\) qui apprécie si la prise en compte de alti et I(alti^2) a permis une hausse significative du \(R^2\) par rapport au modèle gaussien.

Visualisation des prédictions

mu <- modRegQuad_Querc_BC$coefficients[1]; bet <- modRegQuad_Querc_BC$coefficients[2]
gam <- modRegQuad_Querc_BC$coefficients[3]
sigRes <- sd(modRegQuad_Querc_BC$residuals)
plot(modRegQuad_Querc_BC$model[,1]~modRegQuad_Querc_BC$model$alti)
lines(seq(0,1000,1),mu+bet*seq(0,1000,1)+gam*seq(0,1000,1)^2,lwd=2,col=4)

On peut revenir dans l’échelle naturelle :

plot(DBH~alti,data=dataIni_Querc)
lines(seq(0,1000,1),10^(mu+bet*seq(0,1000,1)+gam*seq(0,1000,1)^2),lwd=2,col=2)
lines(seq(0,1000,1),10^(mu+bet*seq(0,1000,1)+gam*seq(0,1000,1)^2+1.96*sigRes),lwd=1,col=2,lty="dashed")
lines(seq(0,1000,1),10^(mu+bet*seq(0,1000,1)+gam*seq(0,1000,1)^2-1.96*sigRes),lwd=1,col=2,lty="dashed")

Tests sur modèles linéaires gaussiens emboîtés

Pour terminer on peut tenter de répondre à la question qui a motivé la section sur le modèle de régression linéaire : est-ce qu’en prenant en compte l’altitude des arbres, on a finalement une explication suffisante de leurs diamètres, où y a-t-il besoin de prendre également en compte le relevé comme on l’avait fait dans la section sur le modèle ANOVA ?

Pour évaluer ce point, on remettre à profit le \(F\)-test vu précédemment. En effet, le \(F\)-test permet également de tester si l’ajout de covariables supplémentaires dans un modèle améliore le coefficient de détermination \(R^2\) de façon significativement plus grande qu’un simple effet d’échantillonnage. On parle de test de modèle “emboîté” car on compare un modèle simple avec un modèle étendu dans lequel on a ajouté des covariables en plus. L’hypothèse nulle est que les covariables ajoutées en plus dans le modèle étendu ne servent à rien, et que l’augmentation du \(R^2\) n’est qu’un pur effet d’échantillonnage.

Dans notre exemple, le test se met en oeuvre comme suit. On commence par définir un super-modèles dans lequel on explique les moyennes des sous-populations à la fois par l’altitude et le relevé. Formellement, ce super modèle s’écrit :

\[ Y_{ij}=\mu_j + \beta_1 \text{alti}_{ij}+ \beta_2 \text{alti}_{ij}^2 +E_{ij} \]

Il est modélisé dans \(R\) comme suit:

modAnovRegQuad_Querc_BC <- lm(I(log10(DBH))~alti+I(alti^2)+releve,data=dataIni_Querc) 

Noter qu’on a réutilisé la transformation Box-Cox retenue dans le modèle de régression linéaire pour bien avoir des modèles emboîtés.

Ensuite on met en oeuvre le \(F\)-test de comparaison de ces deux modèles à l’aide de la fonction anova de R:

anova(modRegQuad_Querc_BC,modAnovRegQuad_Querc_BC)

## Analysis of Variance Table
## 
## Model 1: I(log10(DBH)) ~ alti + I(alti^2)
## Model 2: I(log10(DBH)) ~ alti + I(alti^2) + releve
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1   1458 26.8427                                  
## 2   1450  9.8151  8    17.028 314.44 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

On peut lire que la statistique \(F\) calculée vaut \(106.07\) et que l’augmentation de \(R^2\) associée est très significative (p-valeur inférieure à \(2.2e^{-16}\)). On rejette donc \(H_0\) : l’ajout de la variable relevé augmente encore beaucoup le \(R^2\), même après avoir pris en compte l’altitude et son effet quadratique.

A l’inverse, on peut se demander si une fois qu’on a pris en compte le relevé, ajouter l’altitude et son effet quadratique apporte quelque chose en plus. Pour ce faire on peut tenter de rajouter l’altitude dans le modèle ANOVA développé avec le relevé plus haut. ON reconsidère un super modèle mais avec la variable transformée du modèle ANOVA présenté plus haut

modAnovRegQuad2_Querc_BC <- lm(I(log10(DBH))~alti+I(alti^2)+releve,data=dataIni_Querc) 

et on veut le comparer avec le modèle avec la variable relevé seul, qui est le modèle ANOVA développé plus haut. On doit d’abord réajuster le modèle ANOVA et le super modèle sur le même jeu de données car l’altitude n’est pas renseignée sur certaines observations :

dataIni_Querc2 <- dataIni_Querc[!is.na(dataIni_Querc$alti),]
modAnovDif2_Querc_BC <- lm(I(log10(DBH))~releve,data=dataIni_Querc2)
modAnovRegQuad2_Querc_BC <- lm(I(log10(DBH))~alti+I(alti^2)+releve,data=dataIni_Querc2) 

Ensuite on peut comparer les deux modèles :

anova(modAnovDif2_Querc_BC,modAnovRegQuad2_Querc_BC)

## Analysis of Variance Table
## 
## Model 1: I(log10(DBH)) ~ releve
## Model 2: I(log10(DBH)) ~ alti + I(alti^2) + releve
##   Res.Df    RSS Df Sum of Sq    F   Pr(>F)   
## 1   1452 9.9003                              
## 2   1450 9.8151  2  0.085154 6.29 0.001906 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Le \(F\)-test nous dit que la valeur de la statistique \(F\) vaut \(4.1104\) et qu’il s’agit d’une augmentation significative (\(p=0.01656<0.05\)). Il y a donc un intérêt à ajouter l’altitude et son effet quadratique en plus du relevé. Il faudrait néanmoins vérifier que les hypothèses sur les erreurs sont bien vérifiées dans ce super modèle.

Modèles linéaires mixtes

Motivation pour introduire des modèles mixtes

Dans plusieurs des modèles précédents, la prise en compte d’un effet ‘relevé’ sur la moyenne des sous-populations a souvent donné de bons résultats en termes de vérification des hypothèses et de coefficient de détermination. Cependant, ce choix présente des limites. Notamment il nous empêche de comprendre les facteurs qui peuvent générer ces différences de moyennes entre relevés. Par exemple, imaginons que l’on pense que les différences de latitude entre relevés jouent un rôle dans ces différences. On fabrique une variable de latitude standardisée du relevé que l’on ajoute au jeu de données :

vLatRel <- scale(coef(lm(lati~0+releve,data=dataIni_Querc)))
rownames(vLatRel) <- sub("releve","",rownames(vLatRel))
dataIni_Querc$latiRel <- vLatRel[
  match(dataIni_Querc$releve,rownames(vLatRel)),1]

Chaque relevé a une latitude standardisée :

par(mfrow=c(1,2))
plot(lati~long, data=dataIni,xlab="Longitude WGS84",ylab="Latitude WGS84",pch=3)
text(coorPla[1,],coorPla[2,],colnames(coorPla),col=4,cex=0.5,font=2)
barplot(t(vLatRel)[,order(t(vLatRel),decreasing=TRUE)],las=2,ylab="Latitude standardisée")

Les arbres d’un même relevé sont associés à la même valeur de latitude du relevé :

plot(sort(dataIni_Querc$latiRel),xlab="Arbres",ylab="Latitude standardisée du relevé")

Pour voir si l’effet du relevé est induit en partie par la latitude du relevé, on pourrait penser intègrer la latitude du relevé dans le modèle précédent :

modAnovRegQuadLatRel_Querc_BC <- lm(
  I(log10(DBH))~alti+I(alti^2)+releve+latiRel,data=dataIni_Querc)
summary(modAnovRegQuadLatRel_Querc_BC)

## 
## Call:
## lm(formula = I(log10(DBH)) ~ alti + I(alti^2) + releve + latiRel, 
##     data = dataIni_Querc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.41746 -0.04990 -0.00223  0.05156  0.33300 
## 
## Coefficients: (1 not defined because of singularities)
##                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   2.289e+00  2.733e-01   8.376  < 2e-16 ***
## alti         -1.845e-03  1.452e-03  -1.271 0.204099    
## I(alti^2)     1.065e-06  1.967e-06   0.541 0.588349    
## releveBLO_12 -4.182e-01  4.792e-02  -8.726  < 2e-16 ***
## releveBLO_13  1.170e-01  1.330e-02   8.803  < 2e-16 ***
## releveBLO_17 -4.251e-01  5.331e-02  -7.975 3.06e-15 ***
## releveBLO_21 -2.512e-01  7.274e-02  -3.454 0.000569 ***
## releveBLO_24 -1.651e-01  5.739e-02  -2.877 0.004078 ** 
## releveBLO_27 -2.133e-01  6.617e-02  -3.223 0.001296 ** 
## releveBLO_4  -1.842e-01  4.796e-02  -3.841 0.000128 ***
## releveBLO_9  -1.572e-01  5.187e-02  -3.031 0.002478 ** 
## latiRel              NA         NA      NA       NA    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.08227 on 1450 degrees of freedom
##   (6 observations effacées parce que manquantes)
## Multiple R-squared:  0.7513, Adjusted R-squared:  0.7496 
## F-statistic: 438.1 on 10 and 1450 DF,  p-value: < 2.2e-16

L’effet de la latitude du relevé n’a pas pu être estimé. Cela vient du fait que l’effet ‘relevé’ dans le modèle est déjà par définition suffisant pour absorber toutes les variations entre les sous-populations. Il ne reste aucune variation entre relevés que l’on puisse tenter d’expliquer à l’aide de la latitude.

On peut retirer l’effet relevé du modèle, et ne garder que l’effet latitude :

modRegQuadLatRel_Querc_BC <- lm(
  I(log10(DBH))~alti+I(alti^2)+latiRel,data=dataIni_Querc)
summary(modRegQuadLatRel_Querc_BC)

## 
## Call:
## lm(formula = I(log10(DBH)) ~ alti + I(alti^2) + latiRel, data = dataIni_Querc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.43633 -0.08916 -0.00955  0.07975  0.48120 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  4.022e+00  8.446e-02   47.62   <2e-16 ***
## alti        -1.396e-02  4.681e-04  -29.83   <2e-16 ***
## I(alti^2)    1.879e-05  6.326e-07   29.70   <2e-16 ***
## latiRel      7.725e-02  4.917e-03   15.71   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1255 on 1457 degrees of freedom
##   (6 observations effacées parce que manquantes)
## Multiple R-squared:  0.4184, Adjusted R-squared:  0.4172 
## F-statistic: 349.4 on 3 and 1457 DF,  p-value: < 2.2e-16

On examine l’ajustement de ce modèle :

mod <- modRegQuadLatRel_Querc_BC

par(mfrow=c(1,3))
#DISTRIBUTION SYMETRIQUE
qqPlot(mod,distribution="norm",line="none")

## 354 987 
## 104 419

#HOMOSCEDASTICITE
plot(mod,which=3,pch=3, add.smooth = FALSE)
#Ligne horizontale attendue
abline(h=0.8,col=4,lwd=2)
#Moyenne glissante
lo <- loess(
  sqrt(abs(rstandard(mod)))~mod$fitted.values) 
vFit <- sort(unique(mod$fitted.values))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
#INDEPENDANCE
x <- mod$model$latiRel
plot(mod$residuals~x)
abline(h=0,col=4)
#Moyenne glissante
lo <- loess(
 mod$residuals~x) 
vFit <- sort(unique(x))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")

On constate que les hypothèses du modèles ne sont pas vérifiées. Le troisième diagramme nous montre qu’il y a un effet du relevé qui perdure dans les résidus, qui n’est pas expliqué par la latitude. Dans ce modèle, on a forcé l’effet ‘relevé’ à n’être lié qu’à la latitude du relevé et rien d’autre, ce qui ne semble pas être le cas.

On souhaiterait une solution intermédiaire : expliquer une partie de l’effet du relevé par la latitude tout en autorisant une fluctuation résiduelle de l’effet du relevé non prédite par la latitude. Pour ce faire, on peut proposer une modélisation statistique de l’effet du relevé, c’est à dire considérer que les moyennes des sous populations d’arbres de chaque relevé sont elles-même issues d’une population statistique

Pour un arbre \(i\) dans un relevé \(j\), le dernier modèle linéaire de la section précédente s’écrivait : \[ DBH_{ij}=\mu_j + \beta_1 \text{alti}_{ij}+ \beta_2 \text{alti}_{ij}^2 +E_{ij} \] avec \(E_{ij}\sim\mathcal{N}(0,\sigma^2)\). Dans ce modèle, le paramètre \(\mu_j\) est la moyenne de la sous-population associée au relevé \(j\). Ce dont on parle ici, c’est de considérer que le relevé \(j\) a été échantillonné dans une population de relevés. Une hypothèse simple sur cette population qui intègre un effet de la latitude est de supposer que : \[ \mu_j = \mu+\gamma\,\text{lati}_j+Z_j \] avec \(Z_j\sim\mathcal{N}(0,\tau^2)\)

En d’autres termes, un relevé \(j\) à une latitude \(\text{lati}_j\) est issu d’une population de relevés dont les diamètres moyens sont distribués aléatoirement autour d’une valeur de \(\mu+\gamma\,\text{lati}_j\) selon une gaussienne de variance \(\tau^2\). On reconnaît là un modèle semblable à celui de la régression linéaire, mais appliqué sur l’effet relevé, au lieu de l’être sur les diamètres observés.

Le modèle du diamètre des arbres devient alors: \[ DBH_{ij}=\mu+\gamma\,\text{lati}_j+Z_j + \beta_1 \text{alti}_{ij}+ \beta_2 \text{alti}_{ij}^2 +E_{ij} \] où on a simplement injecté le modèle de \(\mu_j\) dans le modèle de \(DBH_{ij}\). Dans la somme qui constitue le membre de droite de ce modèle, il y a des termes qui sont déterministes : \(\mu\), \(\gamma\,\text{lati}_j\), \(\beta_1 \text{alti}_{ij}\), \(\beta_2 \text{alti}_{ij}^2\). Par déterministe, on entend qu’ils sont complètement connus une fois fixés les paramètres et les observations. On appelle ces termes les ‘effets fixes’ du modèle. On identifie également dans la somme le termes \(E_{ij}\) qui était déjà présent dans les modèles linéaires précédents, et qui est l’erreur associée au diamètre d’un arbre. La nouveauté de cette section est le terme \(Z_j\), qui fait partie de l’effet du relevé, mais est aléatoire (sa valeur n’est pas connue une fois fixés les paramètres et les observations). On parle d’un ‘effet aléatoire’ du relevé. Ici, cet effet aléatoire correspond à la part de l’effet du relevé qui n’est pas associée à la latitude. Ce type de modèle qui mélange des effets fixes et des effets aléatoires est appelé ‘modèle mixte’.

On ajuste ce modèle dans \(R\) de la façon suivante :

library(lme4) #Une bibliothèque pour ajuster des modèles mixtes

## Le chargement a nécessité le package : Matrix

modMixRel <- lmer(I(log10(DBH))~1+alti+I(alti^2)+latiRel+(1|releve),data=dataIni_Querc)

## Warning: Some predictor variables are on very different scales: consider
## rescaling

On a un avertissement car la bibliothèque en question demande que les variables explicatives quantitatives soient standardisées. On crée des variables additionnelles standardisées et on réajuste le modèle :

dataIni_Querc <- data.frame(dataIni_Querc,
  alti.s=scale(dataIni_Querc$alti),
  alti.s2=scale(dataIni_Querc$alti)^2)
modMixRel <- lmer(I(log10(DBH))~1+alti.s+alti.s2+latiRel+(1|releve),data=dataIni_Querc)

Dans la formule \(R\), l’effet aléatoire \(Z_j\) est codé par ‘(1|releve)’.

Diagnostic des hypothèses

Comme précédemment, la première étape de travail sur le modèle est la vérification de l’adéquation des hypothèses du modèle aux données. On a deux séries d’hypothèse dans les modèle mixtes : les hypothèses sur la distribution des effets aléatoire (\(Z_j\)) et les hypothèses sur les erreurs des observations (\(E_{ij}\)).

Diagnostics sur les erreurs résiduelles : on fait comme avant

Les hypothèses sur les erreurs résiduelles \(E_{ij}\) sont les mémes que pour les sections précédentes : normalité, homoscédasticité et indépendance.

Cependant, pour obtenir des résidus associés à un arbre \(i\) dans un relevé \(j\), il faudrait connaître l’effet aléatoire \(Z_j\) associé au relevé \(j\). Cet effet est aléatoire, on ne le connaît pas et on ne l’estime pas directement au moment de l’ajustement du modèle (on estime la variance de la population des effets aléatoires \(\tau^2\)).

IL faut donc construire un estimateur de \(Z_j\) a posteriori, c’est à dire dans un second temps, une fois l’ajustement du modèle effectué. Les outils usuels d’ajustement des modèles mixtes disposent de fonctions pour effectuer des estimations a posteriori des \(Z_j\). On note \(\hat{Z}_j\) l’estimateur de \(Z_j\). Dans lme4, on obtient les estimateurs \(\hat{Z}_j\) comme suit :

ranef(modMixRel)

## $releve
##        (Intercept)
## BLO_1   0.08480769
## BLO_12 -0.21701308
## BLO_13  0.29879887
## BLO_17 -0.19932093
## BLO_21 -0.01967755
## BLO_24  0.07713684
## BLO_27  0.04201704
## BLO_4  -0.05585907
## BLO_9  -0.01088983
## 
## with conditional variances for "releve"

Grâce à ces estimateurs a posteriori des \(\hat{Z}_j\), on peut obtenir des résidus \(\hat{E}_{ij}\) pour chaque arbre : \[ \hat{E}_{ij}=Y_{ij}-( \hat{\mu}+\gamma \, \text{lati}_j+\hat{Z}_j +\hat{\beta}_1\text{alti}_{ij} +\hat{\beta}_2\text{alti}_{ij}^2 ) \]

On appelle ces résidus des ‘résidus conditionnels’, car ils sont calculés conditionnellement à la valeur estimée des effets aléatoires. Dans lme4, on obtient les résidus conditionnels comme suit :

residuals(modMixRel)

On peut effectuer les trois diagnostics usuels sur ces résidus :

par(mfrow=c(1,3))
#DISTRIBUTION SYMETRIQUE
qqPlot(residuals(modMixRel),distribution="norm",line="none")

## 370 987 
## 117 419

#HOMOSCEDASTICITE
vAmp <- sqrt(abs(scale(residuals(modMixRel))))
plot(predict(modMixRel),vAmp)
abline(h=0.8,col=4,lwd=2)
#Moyenne glissante
lo <- loess(vAmp~predict(modMixRel)) 
vFit <- sort(unique(predict(modMixRel)))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")

#INDEPENDANCE
plot(residuals(modMixRel)~predict(modMixRel))
abline(h=0,col=4)
#Moyenne glissante
lo <- loess(residuals(modMixRel)~predict(modMixRel)) 
vFit <- sort(unique(predict(modMixRel)))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")

On ne constate rien d’alarmant sur ces diagnostics.

Diagnostic sur les effets aléatoires

Il reste à évaluer les hypothèses du modèle sur les effets aléatoires \(Z_j\). Il s’agit d’hypothèses très similaires à celles faites sur les erreurs des observations : normalité, homoscédasticité et indépendance. On effectue ces diagnostics directement sur les \(\hat{Z}_j\). On commence par le diagnostic de normalité :

qqPlot(ranef(modMixRel)$releve[,1],distribution="norm",line="none")

## [1] 3 2

Pas de problème marqué, mais on a peu de points pour juger. On regarde l’homoscédasticité en fonctione de la latitude du relevé :

#HOMOSCEDASTICITE
vAmp <- sqrt(abs(scale(ranef(modMixRel)$releve[,1])))
vX <- modMixRel@frame$latiRel[
    match(rownames(ranef(modMixRel)$releve),modMixRel@frame$releve)]
plot(vX,vAmp,ylim=c(0,2))
abline(h=0.8,col=4,lwd=2)
#Moyenne glissante
lo <- loess(vAmp~vX) 
vFit <- sort(unique(vX))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")

Ici non plus, pas de problème à noter, mais on a peu de points pour en juger, donc peu de puissance pour détecter d’éventuelles anomalies. On a regardé l’homoscédasticité le long de la latitude du relevé. On aurait pu faire d’autres choix, mais celui-ci se rapproche le plus des diagnostics en fonction de la valeur prédite dont on a l’habitude dans les modèles linéaires classiques.

Enfin on regarde l’indépendance, ici encore en fonction de la latitude du relevé :

#INDEPENDANCE
plot(ranef(modMixRel)$releve[,1]~vX,ylim=c(-2,2))
abline(h=0,col=4)
#Moyenne glissante
lo <- loess(ranef(modMixRel)$releve[,1]~vX) 
vFit <- sort(unique(vX))
predLo <- predict(lo,vFit,se=TRUE)
lines(predLo$fit~vFit,col=2,lwd=2)
#Enveloppe de confiance autour de cette moyenne glissante :
nFit <- length(vFit); ICBonf <- qnorm(1-0.05/2/nFit)
lines(predLo$fit+ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")
lines(predLo$fit-ICBonf*predLo$se.fit~vFit,col=2,lwd=2,lty="dashed")

On fait les mêmes constats, pas de problème détecté, mais peu de puissance.

Analyse des sorties

Le résumé de l’ajustement du modèle s’obtient via la fonction usuelle :

summary(modMixRel)

## Linear mixed model fit by REML ['lmerMod']
## Formula: I(log10(DBH)) ~ 1 + alti.s + alti.s2 + latiRel + (1 | releve)
##    Data: dataIni_Querc
## 
## REML criterion at convergence: -3076.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -5.0606 -0.6075 -0.0155  0.6367  4.0366 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  releve   (Intercept) 0.02844  0.16864 
##  Residual             0.00677  0.08228 
## Number of obs: 1461, groups:  releve, 9
## 
## Fixed effects:
##              Estimate Std. Error t value
## (Intercept)  1.605693   0.057302  28.022
## alti.s      -0.063061   0.019756  -3.192
## alti.s2      0.006221   0.007474   0.832
## latiRel      0.049926   0.061073   0.817
## 
## Correlation of Fixed Effects:
##         (Intr) alti.s alt.s2
## alti.s   0.038              
## alti.s2 -0.184 -0.412       
## latiRel  0.001 -0.206  0.040

La table de résumé sur les effets aléatoires nous indique que la variance (\(\tau^2\)) de l’effet aléatoire \(Z_j\) du relevé est estimée à 0.02844. La variance de l’erreur associée aux observations est quant à elle estimée à 0.00677.

La table sur les effets fixes nous donne les estimateurs de l’ordonnée à l’origine (Intercept, \(\hat{\mu}\)), le pente estimée de l’effet de l’altitude (alti.s, \(\hat{\beta_1}\)), le pente estimée de l’effet de l’altitude au carré (alti.s2, \(\hat{\beta_2}\)), la pente estimée de l’effet de la latitude du relevé sur l’effet relevé (latiRel, \(\gamma\)). Pour chaque estimateur, on a une erreur standard qui permet de déterminer des intervalles de confiance pour les ‘vraies’ valeurs de paramètres.

Les intervalles de confiances s’obtiennent toujours avec la fonction :

confint(modMixRel)

## Computing profile confidence intervals ...

##                    2.5 %      97.5 %
## .sig01       0.093088878  0.25579799
## .sigma       0.079331857  0.08532121
## (Intercept)  1.497521511  1.71504384
## alti.s      -0.100101134 -0.02192809
## alti.s2     -0.007897805  0.02167766
## latiRel     -0.065161325  0.16664898

On note que dans le cadre du package lme4, la fonction fournit des intervalles de confiances des effets fixes mais également des intervalles de confince des écarts-types des effets aléatoires (\(\tau\), ligne .sig01 ci-dessus) et erreurs associées aux observations (\(\sigma\), ligne .sigma ci-dessus)

On peut regarder l’effet qu’a eu l’introduction de la modélisation de l’effet relevé en fonction de la latitude sur l’estimation des autres effets fixes. Pour ce faire, on compare les estimateurs des autres effets fixes (alti.s et alti.s2) entre modèle mixte et modèle avec relevé en effet fixe. On commence par réajuster le modèle avec relevé en effet fixe sur exactement les mêmes données que le modèle mixte, avec les altitudes standardisées :

modAnovRegQuad_Querc_BC_dataMix <- lm(`I(log10(DBH))`~0+alti.s+alti.s2+releve,data=modMixRel@frame)

print("Modèle mixte :")

## [1] "Modèle mixte :"

confint(modMixRel)[c("alti.s","alti.s2"),]

## Computing profile confidence intervals ...

##                2.5 %      97.5 %
## alti.s  -0.100101134 -0.02192809
## alti.s2 -0.007897805  0.02167766

print("Modèle avec relevé en effet fixe :")

## [1] "Modèle avec relevé en effet fixe :"

confint(modAnovRegQuad_Querc_BC_dataMix)[c("alti.s","alti.s2"),]

##               2.5 %      97.5 %
## alti.s  -0.11181802 -0.03063124
## alti.s2 -0.01078826  0.01901248

On constate qu’on a des intervalles de confiance très similaires entre les deux modèles, ce qui est rassurant : la modélisation de l’effet relevé introduite dans le modèle mixte n’a pas perturbée la perception de l’effet altitude sur les arbres. Notamment on a bien dans les deux cas que l’effet de alti.s est négatif à un niveau de confiance de 95%, tandis que l’intervalle de confiance à 95% de alti.s2 recouvre 0, suggérant que cette pente n’est pas significativement différent de 0.

Contrairement aux fonctions vues plus haut pour les modèles linéaires, le résumé du modèle fournit par lme4 ne donne pas de p-valeurs associées au fait de tester si chaque effet fixe du modèle est différent de 0. Néanmoins on peut faire ce test à la main. Pour un effet en particulier, on peut utiliser un test de modèle emboîté via la fonction anova vue précédemment. Par exemple, on pressent que alti.s2 n’a pas d’effet significatif sur le diamètre des arbres, on le teste comme suit :

modMixRel_noAlti2 <- lmer(I(log10(DBH))~1+alti.s+latiRel+(1|releve),data=dataIni_Querc)
anova(modMixRel_noAlti2,modMixRel,test="Chisq")

## refitting model(s) with ML (instead of REML)

## Data: dataIni_Querc
## Models:
## modMixRel_noAlti2: I(log10(DBH)) ~ 1 + alti.s + latiRel + (1 | releve)
## modMixRel: I(log10(DBH)) ~ 1 + alti.s + alti.s2 + latiRel + (1 | releve)
##                   npar     AIC     BIC logLik deviance  Chisq Df Pr(>Chisq)
## modMixRel_noAlti2    5 -3088.1 -3061.7 1549.0  -3098.1                     
## modMixRel            6 -3086.9 -3055.2 1549.5  -3098.9 0.8387  1     0.3598

On obtient une p-valeur de 0.36, ce qui nous conforte dans le caractère non significatif de l’effet de alti.s2. On peut plus systématiquement tester tous les effets fixes du modèle mixte via la fonction :

drop1(modMixRel,test="Chisq")

## Single term deletions
## 
## Model:
## I(log10(DBH)) ~ 1 + alti.s + alti.s2 + latiRel + (1 | releve)
##         npar     AIC    LRT Pr(Chi)   
## <none>       -3086.9                  
## alti.s     1 -3079.4 9.4955 0.00206 **
## alti.s2    1 -3088.1 0.8387 0.35976   
## latiRel    1 -3088.1 0.8068 0.36906   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Cette opération recrée en quelque sorte la colonne de p-valeurs manquante dans le résumé du modèle. On constate notamment que l’effet de la latitude sur le relévé n’est pas significatif : un modèle gaussien de l’effet relevé sans covariable aurait tout aussi bien convenu.

Pour terminer, on peut souhaiter obtenir une quantification de la capacité globale du modèle à expliquer les variations de diamètre entre les arbres. Dans les modèles classiques, le coefficient de détermination \(R^2\) joue ce rôle, en quantifiant la part de variance expliquée par le modèle, ou autrement dit la part d’erreur quadratique absobrée par le modèle par rapport à un modèle sans covariables ni sous-populations. La notion de coefficient de détermination dans un modèle mixte est ambigüe. En effet, il faut décider comment on considère les effets aléatoires \(Z_j\) : est-ce qu’on les compte dans l’absorption de l’erreur ou non ? Les deux options sont justifiables. D’un côté, on dispose d’estimateurs a posteriori de ces effets (les \(\hat{Z}_j\)), ce qui tendrait à dire qu’on les a quantifiés et qu’on peut les compter dans la partie explicative du modèle. D’un autre côté, ces estimations n’émanent pas directement du modèle (qui n’estime que \(\tau^2\)) mais sont faites dans un second temps, en utilisant une deuxième fois les données. On peut donc discuter le fait d’attribuer leur pouvoir explicatif au modèle lui même. On a donc logiquement deux valeurs possibles du \(R^2\) dans un modèle mixte. Le \(R^2\) conditionnel considère que les effets aléatoires contribuent à l’explication, tandis que le \(R^2\) marginal considère qu’ils n’y contribuent pas. Logiquement, le \(R^2\) marginal est toujours plus faible que le \(R^2\) conditionnel.

Pour obtenir ces valeurs de \(R^2\), on utilise le package MuMIn :

library(MuMIn)
r.squaredGLMM(modMixRel)

## Warning: 'r.squaredGLMM' now calculates a revised statistic. See the help page.

##             R2m       R2c
## [1,] 0.05537556 0.8183667

On obtient un \(R^2\) marginal de 0.05 et un \(R^2\) conditionnel de 0.82. La différence est forte, suggérant deux choses. Tout d’abord, les fluctuations inter-relevé des moyennes de sous-populations sont très fortes, ce qu’on savait déjà vu le \(R^2\) du le modèle avec relevé en effet fixe, et vu la forte valeur estimée de la variance de l’effet aléatoire relevé dans le modèle mixte. Ensuite, la latitude du relevé contribue peu à absorber cette effet car la latitude du relevé contribue au \(R^2\) marginal, qui est faible. Ici encore, on pouvait s’y attendre vu le caractère non significatif de la latitude du relevé obtenue lors des tests systématique par le drop1.

Modèles linéaires généralisés

Une limite de l’hypothèse gaussienne

Pour certains modèles, l’hypothèse gaussienne des résidus ne sera pas vérifiée, même après transformation. Par exemple, si on étude une variable binaire décrivant la présence (1) ou l’absence (0) d’une cavité sur la base du tronc de l’arbre, il est peu probable que des résidus montre une forme en cloche. On va vérifier ce point. On fabrique la variable binaire relevant la présence d’une cavité sur le tronc:

dataIni_Querc$cavPA <- (dataIni_Querc$cav_basses_presence_cavites=="oui")

On applique un modèle linéaire en fonction de l’altitude et du relevé et on évalue l’hypothèse de normalité des erreurs :

modCav <- lm(cavPA~alti+releve,data=dataIni_Querc)
qqPlot(modCav,distribution="norm",line="none")

## 129 171 
##  19  58

On observe que le QQ-plot montre une discontinuité marquée et ne suit pas du tout la diagonale pour les valeurs intermédiaires. On peut tenter une transformation. Les transformation Box-Cox ne s’appliquent que sur des variables strictement positives, on translate donc la variable binaire avant de transformer (en ajoutant 1 à toutes les valeurs) :

modCav2 <- lm(I(1+cavPA)~alti+releve,data=dataIni_Querc)
modCav2_pT <- powerTransform(modCav2)
modCav2_BC <- lm(I((1+cavPA)^modCav2_pT$roundlam)~alti+releve,data=dataIni_Querc)
qqPlot(modCav2_BC,distribution="norm",line="none")

## 129 171 
##  19  58

On récupère le même QQ-plot à l’envers (parce que la transformation Box Cox retenue ici est décroissante). On n’a donc rien réglé en transformant, parce que la variable transformée ne prend toujours que deux valeurs, c’est donc en gros toujours la même variable binaire à un facteur d’échelle près qui ne change rien au problème…

On va bien être obligé de changer de modèle…

Le modèle linéaire généralisé binomial

Présentation

Pour des variables binaires comme la présence d’une cavité, il faut renoncer à l’hypothèse gaussienne pour se tourner vers le modèle de Bernoulli. Ce modèle postule que la variable réponse est distribuée suivant une loi de Bernouilli (vaut \(1\) avec une certaine probabilité \(p\) et \(0\) avec une probabilité \(1-p\)) dans la population sous-jacente à l’échantillon. On suppose toujours que les observations sont tirées de façon indépendante dans la population et on note tout cela mathématiquement :

\(\text{PA}_i \sim \mathcal{B}(p)\)

Si on veut poser un modèle où l’altitude des arbres et le relevé définissent des sous populations, on pourrait essayer de mimer ce qu’on a fait dans le cadre du modèle linéaire gaussien, à savoir postuler que la probabilité de présence d’une cavité pour un arbre \(i\) dans un relevé \(j\) avec une altitude \(\text{alti}_{ij}\) s’écrit:

\[p_{ij} = \mu_j+\gamma \,\text{alti}_{ij}\]

Le problème c’est que le membre de droite de cette définition est une droite, qui peut devenir negative ou supérieure à 1 si \(\text{alti}_{ij}\) prend des valeurs trop extrêmes, ce qui donnerait lieu à une erreur car une porbabilité est par définition comprise entre 0 et 1. Pour ce prémunir de ce bug, on transforme les probabilités \(p_{ij}\) par une fonction \(\Phi\) pour qu’elles puissent sortir de \(\left[0,1\right]\) et on modélise ces probabilités transformées :

\[\Phi(p_{ij}) = \mu_j+\gamma \,\text{alti}_{ij}\]

Le modèle ci-dessus est appelé modèle linéaire généralisé binomial.

La fonction \(\Phi\) utilisée pour transformer les probabilités est appelée ‘fonction de lien’. Il y a plein de choix possibles pour la fonction de lien. Un choix classique est la fonction ‘logit’: \[\text{logit}(p) = \log\left(\frac{p}{1-p}\right)\] Quand on choisit cette fonction, on parle de régression ‘logistique’. On fera systématiquement ce choix dans la suite.

L’erreur associée aux observations \(E_{ij}\) dans le modèle linéaire généralisé binomial se définit comme suit:

\[E_{ij} = \text{PA}_{ij}-p_{ij} = \text{PA}_{ij}-\Phi^{-1}\left(\mu_j+\gamma \,\text{alti}_{ij}\right)\]

Contrairement au modèle gausien, cette erreur est ici complètement déterminée quand on connait les effets fixes et les observations, on n’a pas d’équivalent du paramètre \(\sigma^2\) des modèles linéaires gaussiens. La distribution de \(E_{ij}\) est binaire : l’erreur vaut \(1-p_{ij}\) avec une probabilité \(p_{ij}\) et vaut \(-p_{ij}\) avec une probabilité \(1-p_{ij}\). On peut vérifier qu’en moyenne cette erreur est bien nulle. Sa variance vaut \(p_{ij}(1-p_{ij})\). On note que la variance de l’erreur qui varie avec \(p_{ij}\) : on n’a pas l’hypothèse d’homoscédasticité dans le modèle linéaire généralisé binomial.

Ajustement et évaluation des hypothèses

On peut ajuster un modèle linéaire généralisé avec l’outil glm :

modBinCavPA <- glm(cavPA~alti+releve,
  data=dataIni_Querc,
  family=binomial(link = "logit"))

Ici on utilise le jeu de donnée où on a retiré les observations avec des NAs présents dans les variables utiles, pour éviter des erreurs dans les opérations à venir.

Tout comme pour les sections précédentes, la première étape consiste en la vérification des hypothèses du modèle, à savoir que : (i) les erreurs ont bien la distribution précisée ci-dessus, et (ii) sont indépendantes. Pour cela on va examiner les résidus du modèle, comme on l’a fait jusqu’ici. Les résidus sont définis ici comme : \[\hat{E}_{ij}=\text{PA}_{ij}-\hat{p}_{ij}\] où \[\hat{p}_{ij} = \Phi^{-1}(\hat{\mu}_j+\hat{\gamma} \,\text{alti}_{ij})\] est la réponse prédite par le modèle pour l’observation \(i\) dans le relevé \(j\) tands que \(\hat{\mu}_j+\hat{\gamma} \,\text{alti}_{ij}\) est appelé le prédicteur linéaire pour l’observation \(i\) das le relevé \(j\).

Les résidus sont donc également des variables binaires, comme les erreurs qu’ils estiment. \(\hat{E}_{ij}\) vaut \(1-\hat{p}_{ij}\) avec une probabilité \(p_{ij}\) et vaut \(-\hat{p}_{ij}\) avec une probabilité \(1-p_{ij}\). Le résidu de l’arbre \(i\) du relevé \(j\) a une moyenne de \(p_{ij}-\hat{p}_{ij}\) qui est le appelée le biais du modèle pour l’observation en question. Contrairement au modèles linéaires gaussiens, le biais peut être non nuls dans les modèles linéaire généralisés, du fait de l’application d’une fonction de lien non linéaire.

Il est difficile de diagnostiquer par un examen visuel direct des résidus d’un modèle linéaire généralisé binomial si les hypothèses sont vérifiées, du fait qu’ils sont binaires, hétéroscédastiques et potentiellement biaisés même quand les hypothèses sont valides.

Une logique de diagnostic visuel pour les modèles linéaires généralisés est donc de transformer les résidus d’une façon qui devraient les rendre continus, homoscédastiques et non-biaisés si les hypothèses du modèle sont valides. Puis de faire des diagnistics visuels sur ces résidus transformés. La bibliothèque DHARMa permet d’effectuer ces transformations :

library(DHARMa)

## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')

Grâce à cette bibliothèque, on peut transformer les résidus :

residusTransfo <- simulateResiduals(modBinCavPA ,n=1000)

En gros ce que fait cette transformation c’est rendre les résidus continus en y ajoutant un bruit aléatoire, puis évaluer comment le résidu bruité se situe dans une distribution de valeurs de résidus obtenue par simulation de données avec le modèle étudié en déterminant le quantile auquel il correspond. Si les hypothèses du modèle sont valides, on s’attend à ce que le quantile du résidu bruité observé dans la distribution des résidus bruités simulés suive une loi uniforme entre \(0\) et \(1\). On utilise le quantile comme résidu transformé.

On peut contrôler que ces résidus ont bien une distribution uniforme:

testUniformity(residusTransfo)

## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  simulationOutput$scaledResiduals
## D = 0.014347, p-value = 0.9235
## alternative hypothesis: two-sided

On peut tester qu’on a bien une distribution uniforme par relevé :

testCategorical(residusTransfo,modBinCavPA$model$releve)

## $uniformity
## $uniformity$details
## catPred: BLO_1
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.065648, p-value = 0.6493
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: BLO_12
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.054791, p-value = 0.28
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: BLO_13
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.083963, p-value = 0.7427
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: BLO_17
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.033404, p-value = 0.8296
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: BLO_21
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.072205, p-value = 0.4969
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: BLO_24
## 
##  Exact one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.063014, p-value = 0.9125
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: BLO_27
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.069993, p-value = 0.6484
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: BLO_4
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.079531, p-value = 0.4554
## alternative hypothesis: two-sided
## 
## ------------------------------------------------------------ 
## catPred: BLO_9
## 
##  Asymptotic one-sample Kolmogorov-Smirnov test
## 
## data:  dd[x, ]
## D = 0.079892, p-value = 0.2339
## alternative hypothesis: two-sided
## 
## 
## $uniformity$p.value
## [1] 0.6492620 0.2800278 0.7427349 0.8296397 0.4968539 0.9125106 0.6483756
## [8] 0.4553868 0.2338599
## 
## $uniformity$p.value.cor
## [1] 1 1 1 1 1 1 1 1 1
## 
## 
## $homogeneity
## Levene's Test for Homogeneity of Variance (center = median)
##         Df F value Pr(>F)
## group    8  0.9959 0.4372
##       1457

Enfin, on peut vérifier que les résidus sont bien répartis de façon homogène quelque soit la valeur prédite ou la valeur d’une variable (ici l’altitude) :

testQuantiles(residusTransfo,predictor=modBinCavPA$model$alti)

## 
##  Test for location of quantiles via qgam
## 
## data:  simulationOutput
## p-value = 0.4625
## alternative hypothesis: both

Ces diagnostics sont une sorte d’équivalent des tests que l’on a mené sur les modèles linéaires.

Analyse des sorties

On analyse les sorties d’un modèle linéaire généralisé d’une façon très similaire à ce qu’on fait pour les modèles linéaires. Le résumé d’obtient comme suit :

summary(modBinCavPA)

## 
## Call:
## glm(formula = cavPA ~ alti + releve, family = binomial(link = "logit"), 
##     data = dataIni_Querc)
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   7.30385    5.12004   1.427 0.153718    
## alti         -0.01742    0.01091  -1.597 0.110253    
## releveBLO_12 -5.28064    1.52047  -3.473 0.000515 ***
## releveBLO_13  0.32300    0.33102   0.976 0.329177    
## releveBLO_17 -6.36139    1.81888  -3.497 0.000470 ***
## releveBLO_21 -4.24114    2.55226  -1.662 0.096569 .  
## releveBLO_24 -3.16702    1.93276  -1.639 0.101296    
## releveBLO_27 -4.45740    2.33086  -1.912 0.055832 .  
## releveBLO_4  -4.18095    1.51093  -2.767 0.005655 ** 
## releveBLO_9  -3.61230    1.68352  -2.146 0.031898 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1068.15  on 1465  degrees of freedom
## Residual deviance:  866.57  on 1456  degrees of freedom
##   (1 observation effacée parce que manquante)
## AIC: 886.57
## 
## Number of Fisher Scoring iterations: 7

On y retrouve les valeurs des effets estimés, les erreurs standards associées, qui permettent de déterminer des intervalles de confiance des ‘vraies’ valeurs des paramètres. On a également un test pour déterminer si ces valeurs dévient significativement de \(0\). A noter que la statistique utilisée pour le test (\(z\)) est différente de celle utilisée dans les modèles linéaires classiques (\(t\)) mais on n’approfondira pas ce point. La dernière colonne donne la p-valeur du test.

Les intervalles de confiance des effets fixes s’obtiennent toujours via :

confint(modBinCavPA)

## Attente de la réalisation du profilage...

##                    2.5 %       97.5 %
## (Intercept)  -2.73483269 17.381087771
## alti         -0.03891814  0.003931876
## releveBLO_12 -8.27399238 -2.300001976
## releveBLO_13 -0.33062299  0.970927118
## releveBLO_17 -9.95230990 -2.804348198
## releveBLO_21 -9.25998574  0.767947446
## releveBLO_24 -6.96405375  0.629857575
## releveBLO_27 -9.04335441  0.114706622
## releveBLO_4  -7.16048380 -1.223893633
## releveBLO_9  -6.92818509 -0.313155188

On note que l’intervalle de confiance de l’altitude coupe \(0\), en conformité avec le test qui évalue si l’effet diffère de \(0\) dans le résumé, qui donne une p-valeur de \(0.11\), supérieure à \(0.05\) donc.

On pourrait souhaiter évaluer si le relevé contribue significativement à expliquer la présence de cavité. POur ce faire, on peut faire un test de modèle emboîté, d’une façon très similaire à celle utilisée dans le cadre des modèles linéaires:

anova(
  glm(cavPA~alti,data=dataIni_Querc,family=binomial(link="logit")),
  modBinCavPA,test="Chisq")

## Analysis of Deviance Table
## 
## Model 1: cavPA ~ alti
## Model 2: cavPA ~ alti + releve
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1      1464    1052.26                          
## 2      1456     866.57  8   185.69 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

A noter qu’il s’agit désormais d’un test avec une statistique du \(\chi^2\) (‘Chisq’) et non plus un test \(F\), mais on interprète de la même façon. Ici on voit que le relevé contribue très fortement à expliquer la présence de cavité.

Reste à donner une quantification globale du potentiel explicatif du modèle. Il y a plusieurs traduction possible du coefficient de détermination dans le cadre des modèles linéaires généralisés. On en donne un ici, le \(R^2\) de McFadden, qui quantifie de combien on a progressé vers le meilleur modèle possible en terme de vraisemblance (probabilité des observations sachant les paramètres) quand on est passé d’un modèle sans covariables au modèle actuel. On calcule ce \(R^2\) à la main:

mod_0 <- glm(cavPA~1,data=dataIni_Querc,family=binomial(link="logit"))
R2_mcFadden <- 1-modBinCavPA$deviance/mod_0$deviance
R2_mcFadden

## [1] 0.1889081